CC BY
44Актуальные проблемы нефти и газа ■ Вып. 2(21) 2018 ■ http://oilgasjournal.ru
ОБРАТНЫЕ ЗАДАЧИ ПО ИДЕНТИФИКАЦИИ ПАРАМЕТРОВ ПЛАСТА
(ЗАДАЧИ HISTORY MATCHING)
Э.С. Закиров, С.Н. Закиров, И.М. Индрупский, Д.П. Аникеев ИПНГ РАН, e-mail: ezakirov@ogri.ru
Современные методы геологии, геофизики, исследований кернов, пластов и скважин позволяют строить трехмерные (3D) геолого-гидродинамические модели продуктивных пластов. Однако, как правило, прогнозные расчеты на основе этих моделей приводят к довольно существенному отличию результатов численного моделирования от фактических данных эксплуатации скважин. Поэтому давно начали проводиться исследования с целью идентификации параметров пласта с использованием фактических показателей разработки [1, 2 и др.].
Уточненная, адаптированная 3D модель пласта позволяет осуществлять достоверный прогноз показателей добычи нефти, газа, воды, конденсата как по отдельным скважинам, так и месторождению в целом на ту или иную перспективу.
В задачах идентификации предполагается, что в течение нескольких лет разработки получены фактические данные об изменении во времени дебитов скважин по нефти, газу и воде, забойных и пластовых давлений по скважинам, насыщенностей фаз вдоль забоев скважин (вскрытых интервалов пласта). Очевидно, что эти данные несут в себе информацию о параметрах разрабатываемой залежи. Поэтому естественно стремление на основе этой информации уточнить емкостные и фильтрационные параметры во всем объеме рассматриваемого пласта.
Соответствующие задачи получили название обратных (коэффициентных) задач, задач идентификации коллекторских свойств (history matching problem).
Обратные задачи формулируются в оптимизационной постановке. В качестве функционала, критерия качества, цели обычно принимается взвешенная сумма среднеквадратичных невязок между расчетными и фактическими значениями всех измеряемых параметров по всем добывающим, нагнетательным и наблюдательным скважинам за рассматриваемый период истории разработки месторождения. В более продвинутых постановках целевой функционал может также включать регуляризирующие слагаемые, в том числе обосновываемые из статистических принципов (например, по методу Байеса).
Алгоритмы решения разных задач идентификации параметров пласта излагаются, например, в работах [3, 4, 5, 6, 7, 8, 9]. Общим для этих алгоритмов является использование градиентной процедуры. Наиболее эффективным подходом к вычислению необходимых производных целевого функционала по искомым параметрам являются методы теории оптимального управления, иначе называемые сопряженными методами. При этом на каждой итерации обратной задачи осуществляется решение двух задач:
• прямой краевой 3Б многофазной задачи с использованием заданных фактических режимов работы скважин по выбранному контрольному показателю, например, дебиту по нефти;
• сопряженной краевой 3Б «многофазной» задачи относительно вспомогательных (сопряженных) функций.
Кроме того, специальная процедура, включающая решение дополнительной краевой задачи, используется для выбора шага смещения вдоль направления поиска. Все это значительно ускоряет целенаправленный поиск минимума целевого функционала вне зависимости от количества уточняемых параметров, числа скважин, объема фактической информации об эксплуатации скважин и продолжительности истории разработки месторождения.
Приведем достаточно общее описание алгоритма автоматизированной адаптации модели к истории разработки, использующего эффективные методы теории оптимального управления для вычисления градиента целевой функции.
Запишем целевую функцию в следующем виде [10, 11]:
Здесь индекс у относится к определенному моменту времени измерения, N - общее число моментов измерения; и у] соответствуют векторам измеренных и вычисленных значений динамических данных (забойные, пластовые давления, дебиты компонентов, соотношения фаз в потоке (обводненность, газовый фактор и т.п.), насыщенности фаз и другие доступные динамические данные в различных скважинах/интервалах/сеточных блоках), соответствующие у'-ому измерению; верхний индекс Т обозначает транспонирование. Симметричная (и часто диагональная) весовая матрица О служит для нормализации и приведения в соответствие доверительных интервалов различных измерений. В байесовском подходе матрица О является обратной матрицей ковариации, а
(1)
критерий (1) модифицируется добавлением дополнительного априорного члена к целевой функции (1), что является тривиальным для учета в рамках рассматриваемого вычислительного алгоритма. Вектор xj обозначает значения фазовых переменных прямой задачи на j-ом временном шаге во всех сеточных блоках. Например, в трехфазном сеточном блоке типичной модели нелетучей нефти к-ът подвектор вектора г' состоит из давления нефтяной фазы, водо- и газонасыщенности сеточного блока модели. Забойные давления вычисляются для скважин, работающих на заданном дебите компонента, и также включаются в вектор xj [12]. Само моделирование скважины подробно рассматривается в статье данного выпуска .
В формуле (1) и - вектор управляющих параметров. Он определяется решением обратной задачи и обычно состоит из значений пористости и проницаемости сеточных блоков модели, или множителей на пористость и проницаемость для данной зоны пласта. Любые другие параметры, например, коэффициенты функциональных зависимостей функций относительных фазовых проницаемостей, локальная геометрия пласта или свойства водоносного горизонта (как сеточного, так и аналитического, например, типа Фетковича) могут также включаться в вектор и .
Прямая задача состоит в решении системы уравнений в частных производных, описывающих фильтрацию нефти, газа и воды в пласте согласно, например, модели нелетучей нефти (black oil), с соответствующими начальными и граничными условиями [13]. После дискретизации систему нелинейных уравнений прямой задачи на j-ом временно шаге можно записать в виде
F3 (х3, ¿Г"1 ,г/) = 0, j = \JV. (2)
Для решения системы (2) используется метод Ньютона, на каждой итерации а которого возникает следующая линейная система уравнений
F1, (хЛа)) дхЛа) = -FJ'(xJ'(a)) , (3)
где Fyj обозначает матрицу Якоби Fj по переменным х', a Sxj(v) - приращение фазовых
переменных на текущей итерации а .
Заметим, что в общем случае временной интервал между двумя последовательными измерениями в критерии (1) может составлять несколько временных
* См. статью Закиров Э.С., Закиров С.Н., Индрупский И.М., Аникеев Д.П. «О представлении скважины в 3D гидродинамической модели» в данном выпуске.
шагов, но в (1) и (2) мы используем один и тот же временной индекс у для простоты изложения.
Обратная задача ставится как задача минимизации целевой функции (1) на допустимой области й е и при условии выполнения уравнений прямой задачи (2). Соответствующая минимизация осуществляется посредством итерационной градиентной процедуры. Квазиньютоновские процедуры типа BFGS или SSVM [14, 15] доказали свою высокую скорость сходимости и надежность для рассматриваемого типа задач. На каждой итерации обратной задачи вектор управляющих параметров корректируется следующим образом:
й<у+1) =й{у)-(3{У)&У), (4)
где б(у) - направление поиска квазиньютоновского метода, а /?(г) - шаг смещения вдоль
направления поиска. Для определения 1У11 требуется вычислить только градиент V./ по управляющим параметрам й. Эффективная процедура расчета V./ строится с применением дискретной формы принципа максимума Понтрягина. Градиент WJ получают на основе решения сопряженной задачи (5) согласно выражению (6):
('•; )> -№+1)гу>]+1 -¿7й]=жI, (5)
Г т
= Х 2 (у1) Пу / V
(6)
]=1
где № - величина /-ого временного шага, а у/' - сопряженный вектор. При этом выполнено следующее «начальное» условие на сопряженный вектор: ^+1=б. Сопряженная система имеет ту же размерность, что и система уравнений прямой задачи, но решается в обратном направлении по времени. При этом она линейна, поэтому время, требуемое для ее решения, в несколько раз меньше, чем время решения нелинейной прямой задачи вне зависимости от числа управляющих параметров. Как только задача (5) решена, а градиент V./ вычислен, направление поиска Г>(у) легко получается согласно пересчетным формулам используемого квазиньютоновского метода
[14, 15].
Обычный способ определения оптимальной величины шага смещения р(У) требует вычисления целевой функции J в ряде точек, т.е. требуется ряд прогонов симулятора (решение серии прямых задач), что является крайне затратным с вычислительной точки
зрения. Специальная приближенная процедура доказала свою высокую эффективность. Она включает решение задачи для вариаций фазовых переменных 8х3
Р1,8х> 7 = (7)
с начальным условием 8х° = б. После решения задачи (7) оптимальная величина шага смещения вычисляется, исходя из линеаризации зависимости у3 (х3 (й),й) и точки минимума квадратичной функции (1), следующим образом:
N т
^(у4х3(й<у)),й<у))- У3 ) 0.ду3ф{уу)
-*- , (В)
^ёузТ0(У>) П 8у3ф
з=1
принимая во внимание зависимость вариаций 8у3 от 8х3. Снова вычислительная стоимость задачи (7) в несколько раз меньше стоимости прогона одной прямой задачи.
Описанный алгоритм реализован в некоммерческом симуляторе нелетучей нефти 81тММеИ®, обладающем дополнительными опциями по автоматизированной адаптации истории разработки, регулированию разработки на основе оптимального распределения расходов по добывающим и нагнетательным скважинам, а также авторского алгоритма решения задач ремасштабирования (upscaling'а) в 3D глобальной нестационарной постановке [16]. В частности, он способен определять проницаемости вдоль трех главных осей тензора проницаемости, включая анизотропию проницаемости вдоль и поперек напластования по данным специализированных ГДИС по 3D и вертикальному гидропрослушиванию [17].
Геометрическая идентификация
Наряду с задачами идентификации значений параметров пласта в ячейках модели (или их множителей по зонам модели), в случае недостаточной изученности геологического строения месторождения, приходится решать два типа задач геометрической идентификации.
Первая задача заключается в идентификации вертикальных отметок структурных поверхностей (или толщин) продуктивного пласта. Решение данной задачи базируется на фактических данных эксплуатации скважин за прошедшие годы разработки, и постановка задачи аналогична предыдущей. Соответствующий алгоритм решения заключается в идентификации отметок глубин и толщин пропластков за счет использования высотных отметок угловых точек ячеек модели в качестве управляющих
параметров й. В работе [8] впервые осуществлена попытка исследования возможности идентификации пластовой геометрии при помощи методов теории оптимального управления, что оправдало себя на тестовых примерах и примере реальной залежи.
Вторая задача связана с идентификацией размеров и конфигурации пластовой водонапорной системы [18, 19]. Здесь также привлекается фактическая информация об эксплуатации скважин и данные контроля за процессом разработки. В роли управляющих параметров и выступают масштабирующие коэффициенты. Они представляют собой множители, позволяющие сжимать или растягивать границы водоносного бассейна вдоль четырех направлений. Как и в иных задачах идентификации, такие «деформации» водоносного бассейна проводятся до достижения минимума целевого функционала.
Примеры практического решения идентификационных задач
В середине 1980-х годов, впервые в России, возможно и в мире, на основе изложенных алгоритмов была создана комплексная, адаптирующаяся, постоянно действующая геолого-математическая модель разработки газового месторождения Медвежье [6]. Здесь задачи прогнозирования разработки и идентификации коллекторских свойств решались в 3D двухфазной (газ-вода) постановке. При этом модель включала в себя и всю систему обустройства промысла.
Алгоритм решения обратных задач в 1990-х годах успешно использован на ряде объектов. В качестве примера приведем некоторые данные применительно к подземному газохранилищу Лаухштад [20, 21] и месторождению Томмелитен Гамма [9]. Более широкому применению алгоритмов решения задач идентификации параметров пласта на многих объектах мешает недостаточное количество и качество замеров показателей эксплуатации скважин.
Подземное газохранилище Лаухштад (Восточная Германия) создано в истощенном газовом месторождении. На рис. 1 приведено объемное изображение продуктивного пласта и его сеточная аппроксимация. Рассматриваемое ПХГ располагает 17 эксплуатационно-нагнетательными скважинами. Надежные замеры текущих пластовых давлений имеются по 5 скважинам. Адаптация геолого-математической модели пласта выполнена для интервала времени в один год.
Рис. 1. Газовое хранилище Lauchstadt
На рис. 2 в качестве примера для одной из скважин дается сопоставление фактических пластовых давлений с рассчитанными на основе исходной геолого-математической модели пласта. Видны значительные расхождения сопоставляемых давлений. На рис. 3 приводятся те же фактические динамики пластового давления, но в сравнении с расчетными динамиками согласно уточненной модели продуктивного пласта. Результаты сопоставления говорят сами за себя.
140
120
5 н а
«Г § 100
п
со а
ггД
80
60
( -9— факт F^^A 1 расч
f
\t
/
200
400
600
800
Время, сут
Рис. 2. Замеренные и расчетные значения пластового давления по скважине Lt053
(начальное приближение)
0
S н а
«Г =
X
<и
¡а а
100
80
60
/ ^Krlujj^
1
1 V
200
400
600
800
Время, сут
Рис. 3. Замеренные и расчетные значения пластового давления по скважине Lt053
(после history matching)
Разработанный алгоритм решения обратных задач был применен для идентификации параметров продуктивного пласта газоконденсатного месторождения Томмелитен Гамма, расположенного в норвежском секторе Северного моря. Объемное изображение этого месторождения приводится на рис. 4. Для рассматриваемой залежи характерно наличие краевой воды. При моделировании данного месторождения была выбрана трехмерная трехфазная (газ, вода, конденсат) модель пласта. Пласт представлен семью продуктивными пропластками, сильно отличающимися по коллекторским свойствам. При этом верхние 4 пропластка отделены слабо проницаемой заглинизированной перемычкой от нижних пропластков. Залежь разрабатывалась в течение 3,5 лет в режиме истощения, что привело к выпадению конденсата в пласте и призабойных зонах скважин. В результате, насыщенность конденсатом достигла от 20 до 30%, что превысило порог начала его подвижности.
На месторождении пробурено 6 эксплуатационных скважин. Все скважины, за исключением скважины 06H, вскрыли все 7 продуктивных пропластков. Скважина 06H вскрыла лишь нижние три пропластка. Наиболее надежной промысловой информацией являлись динамические забойные давления по скважинам, которые и были выбраны в качестве наблюдаемых параметров. На рис. 5 в качестве примера приводятся результаты
0
расчетов забойного давления для скважины 01Н на исходной и адаптированной геолого-математической модели месторождения.
Рис. 4. Объемное изображение и сеточная аппроксимация месторождения Томмелитен Гамма при идентификации параметров пласта
0 400 800 1200
Время, сут
Рис. 5. Динамика забойных давлений для скважины 01H месторождения Томелитен Гамма
На рис. 6 дается объемное изображение и сеточная аппроксимация нефтегазового месторождения, для которого была выполнена идентификация толщин продуктивного пласта в 3D трехфазной постановке задачи. Об успешности реализации алгоритма свидетельствует, например, рис. 7. На этом рисунке приводятся фактические и расчетные зависимости газового фактора от времени, полученные на исходной и уточненной модели месторождения. Близкими оказались также расчетные и замеренные значения забойного давления по скважинам, а также данные об обводненности добываемой продукции.
Рис. 6. Объемное изображение и сеточная аппроксимация исследованной нефтегазовой залежи
220
140
200 250 300 350 400
Время, сут
Рис. 7. Расчетные и фактические динамики газового фактора для скважины OP0
нефтегазового месторождения
Другой пример реализации геометрической идентификации излагается в работах [18, 19]. Здесь выполнен следующий математический эксперимент. В прямоугольном в плане водоносном пласте с непроницаемыми внешними границами разрабатывается залежь 9 вертикальными скважинами. Залежь асимметрично расположена относительно всех границ водоносного пласта. Смоделирован процесс разработки залежи на период 15 лет при заданных дебитах скважин по жидкости. Полученные показатели разработки (пластовые давления, зависимости от времени газового фактора и обводненности добываемой продукции по скважинам) приняты в качестве фактических. С использованием этой информации в 3D трехфазной постановке решена обратная задача идентификации конфигурации и размеров водоносного бассейна, инициированная с измененного начального приближения. Предложенный алгоритм решения задачи через 24 итерации привел к получению данных о размерах водоносного пласта, которые задавались при решении исходной прямой задачи. Результаты превзошли ожидаемую авторами эффективность алгоритма.
Формулировка и алгоритм решения обратной задачи идентификации в геологически согласованной постановке
Процедура адаптации 3D гидродинамической модели нефтегазоносного пласта к истории разработки месторождения является необходимым этапом в методологии 3D компьютерного моделирования. Эта процедура трудоемка и затратна. Хотя имеют место попытки и успехи в комплексировании разнородных, разномасштабных и разнохарактерных данных при построении 3D геологической модели. Ручная или автоматизированная коррекция 3D гидродинамической модели пласта на этапе адаптации часто производится без требуемого учета особенностей построения 3D геологической модели. Как следствие возникает разрыв между методами и способами построения 3D геологической модели и ее искажением при создании 3D гидродинамической модели пласта, а также при адаптации ее к истории разработки.
В развитие методов автоматизированного решения обратных задач идентификации авторами предложен метод согласованной корректировки 3D гидродинамической и 3D геологической модели пласта на основе уточнения параметров геостатистической модели [22, 23, 24]. Например, таких параметров вариограммы, как диапазон, плато и эффект самородка. Получены выражения для соответствующих производных применительно к известным типам аналитических форм вариограмм. Вычисление функциональных производных критерия качества осуществляется при помощи современных методов
теории оптимального управления. При этом обратная задача рассматривается в максимальной общности: с учетом фациальной неоднородности и анизотропии вариограмм. Уточнению также подвергаются параметры корреляционной зависимости «логарифм проницаемости - пористость» для каждой отдельной фации. В конце раздела приводятся синтетические примеры, подтверждающие работоспособность предлагаемого метода решения 3D обратной задачи.
Рассмотрим задачу геологически согласованной адаптации модели. Модифицируем описанную ранее обратную задачу по уточнению значений пористости и проницаемости индивидуальных сеточных блоков за счет записи целевого функционала (1) в следующем виде:
где V - параметры анизотропной вариограммы (диапазон для главного, вспомогательного и вертикального направлений, эффект самородка и азимутальный угол), а а - параметры корреляции «пористость - логарифм проницаемости» (обычно линейная функция с постоянными коэффициентами) для каждой фации.
В отличие от формулировки (1), новая постановка задачи трактуется следующим образом. Распределение значений коэффициентов пористости и проницаемости по ячейкам модели, формирующее вектор й, не является независимым. Оно строится с использованием заданных значений пористости на скважинах (по данным интерпретации геофизических исследований скважин - ГИС) и методов геостатистики, в частности, процедуры кригинга для их распространения на межскважинное пространство [25]. Проницаемость для каждой ячейки рассчитывается по зависимости от пористости. Таким образом, следует решать обратную задачу относительно V и а .
Описанный в предыдущем разделе алгоритм модифицируется для рассматриваемой задачи естественным образом. Единственная трудность, которую необходимо преодолеть, состоит в вычислении производных по V и а вместо и в (6) и (7). Заметим вначале, что дифференцирование любой векторной функции К = Ё(и(у,а)) выполняется по цепному правилу:
Формулировка новой задачи и вычисление производных
(9)
При этом матрица Якоби йэ вычисляется напрямую, поскольку корреляция пористости с проницаемостью выражается явной функцией. Следовательно, вычисление матрицы Якоби й~ представляет главный интерес. Кроме того, следует обсуждать только матричные строки, соответствующие пористостям сеточных блоков, вычисленных кригингом. Вычисление оставшихся строк тривиально, поскольку проницаемости явным образом связаны с пористостями через а. Таким образом, для упрощения формул предположим далее, что и - вектор пористостей сеточных блоков в пределах определенной фации. И применим общее допущение, что функция плотности вероятности для и гауссова (в противном случае необходимо осуществить операцию нормализации).
Следуя хорошо известной процедуре кригинга, пористость определенного
сеточного блока г (г = 1, , где N ь - число сеточных блоков в модели) вычисляется по формуле
N3^ 1=1
где г - порядковый номер статического значения пористости ия., полученного по ГИС, N,1 - общее число статических данных пористости по скважинам, а А - веса кригинга. Поскольку кригинг - точный метод, нет необходимости исключать зафиксированные значения пористости на скважинах из вектора й. Кроме того, условие наилучшей линейной несмещенной оценки приводит к следующему набору линейных уравнений относительно А (система кригинга):
(12)
Здесь уц =/(<х,У,2>я ><х>У>2>я ~ значение вариограммы, вычисленное для
Гц Г1,2 • • Г 1 ' "А ■ ' Г1,г "
Г 2,1 Г 2,2 • 1 А 2 Г 2,г
Г N,1,1 ,2 1 ГNsd ,г
1 1 • 1 0 _ и _ 1
двух точек статических данных с порядковыми номерами / и /2 (/,/2 = 1, ^), а у1г = у(< х,у,2 >я.,< х,у,2 >г,у) - значение вариограммы, вычисленное для /'-ой точки
статических данных (/ = 1, ) и центра сеточного блока г; ¿и - множитель Лагранжа для
_ 31
условия несмещенной оценки = 0.
1=1
Поскольку ив (11) зафиксированы, то производные иг по у линейно связаны соотношением (11) с матрицей Якоби Л~. В свою очередь, неявная связь между Л и у дается системой (12) и используемой моделью вариограммы. Следовательно, неявные производные по каждому управляющему параметру вариограммы у вычисляются из решения следующей системы линейных уравнений:
Ги
Г 2,1
^ ,1
Г1,2 Г2,2
^ ,2
Гц
Г2,1
Л ду, Л
ду,
Л
СУ, дц ду,
дд
СУ,
Е су1
дЕ су1 0
(13)
где ^(Л,у) = 0 - функциональное выражение для /-ой строчки системы (12) (линейное по
Л и нелинейное по у ).
Набор линейных уравнений (13) должен быть решен в каждом сеточном блоке для каждого управляющего параметра вариограммы. Эффективная реализация возможна благодаря особенностям процедуры кригинга. Поскольку количество статических данных умеренное, матрица слева в (12) обращается прямым образом единожды для каждой итерации обратной задачи. Затем (12) и (13) решаются столько раз, сколько требуется, посредством только умножения матрицы на вектор.
Альтернативная формулировка для проницаемостей Представленный вывод формул (9)-(13) применим, если распределение пористости получают кригингом, а проницаемость рассчитывается посредством корреляции пористости и проницаемости. Следует рассмотреть и другую возможность, когда статические значения проницаемости в скважинах берутся из каротажных данных, то есть проницаемость также определяется процедурой кригинга. В этом случае общая процедура решения обратной задачи близка к вышеописанной, за исключением нескольких позиций.
Во-первых, функция распределения проницаемости обычно близка к лог-нормальной, что означает, что вектор й должен включать значения логарифмов проницаемостей. Также вариограммы для пористости и логарифма проницаемости могут быть различными, так что у включает параметры обеих вариограмм для каждой фации.
Кроме того, данные ГИС по проницаемости обычно далеки от истинных пластовых значений эффективной проницаемости. Поэтому целесообразно в обратной задаче ввести множители на проницаемость, применяемые ко всем значениям статической проницаемости в данной скважине. Эти множители включаются в набор управляющих параметров обратной задачи. Таким образом, формулировка (9) все еще справедлива, но а обозначает вектор скважннных множителей проницаемости.
Теперь осталось описать процедуру вычисления производных проницаемостей сеточных блоков относительно множителей скважинных проницаемостей. Значение проницаемости в сеточном блоке г вычисляется процедурой кригинга (11) (где и теперь обозначает логарифм проницаемости). Значения проницаемости в скважинах не фиксируются, а вычисляются в виде произведения их начального значения на значение скважинного множителя проницаемости на текущей итерации обратной задачи. Таким образом:
п„=1п к:=1п(«:)+), (14)
где верхний индекс обозначает номер итерации обратной задачи, а: - множитель проницаемости в скважине н на текущей итерации :, а значение проницаемости к соответствует скважине н. После подстановки (14) в (9) производная кг по ап вычисляется по формуле:
С »¡й \
дк деп д
да да да,.
е'
ХЛп, »й (; Л »й (з Л = е 1 -Х ЛР = К-Х ^Р» , (15)
'=1 V а» ; '=1 V а»
где Р = 1, если значение проницаемости К, соответствует скважине н, и 0 - в противном случае.
Реализация
Описанные процедуры вычисления производных (коэффициентов чувствительности) проницаемостей и пористостей сеточных блоков относительно V и а были реализованы двумя способами: в виде соответствующего модуля в программном продукте 8тММсИ®, и в виде подключаемого расширения для коммерческого пакета 3D геологического моделирования Ре^е(® фирмы 8сИ1итЪе^ег. Данное расширение полезно при вычислении и визуализации 3D карт производных, чтобы помочь специалистам по геологии и разработке лучше понимать влияние параметров вариограмм и корректирующих множителей скважинной проницаемости на распределение пористости и
проницаемости в межскважинном пространстве (в процессе адаптации модели к истории разработки). Например, рис. 8 показывает пример карты производной пористости по диапазону вариограммы в главном направлении, вычисленной на стандартной тестовой модели Рв&в1® месторождения Галфакс по зоне Tarbert-1. Модель содержит 33909 сеточных блоков, 15 скважин и 490 статических скважинных данных для пористости и проницаемости. В данном случае карта представлена для модели с единой фацией, также рассматривалась версия с 3 флювиальными фациями.
Рис. 8. Карта производной пористости по главному диапазону (одна фация, значение диапазона переоценено)
При реализации общего алгоритма решения обратной задачи каждая его часть тщательно тестировалась. Наконец, общая процедура согласованной адаптации истории протестирована на ряде синтетических примеров. В каждом примере прямая задача решалась на заранее заданном временном интервале с «истинными» значениями управляющих параметров для получения «замеренных» динамик давлений и дебитов по скважинам. Затем, после произвольного изменения значений управляющих параметров, процедура решения обратной задачи запускалась с использованием указанных «замеренных» динамических данных. Результаты выполненных разноплановых синтетических тестов доказали, что «истинные» значения каждого отдельного параметра вариограммы и корреляционной петрофизической зависимости могут быть восстановлены
по наблюдаемым динамическим данным. Также на идентифицируемость протестированы пары и тройки параметров. Для рассмотренных тестовых примеров «истинные» значения управляющих параметров были восстановлены с высокой точностью. В качестве примера рис. 9 и 10 демонстрируют динамики «замеренных» забойных давлений для двух скважин, в сравнении с вычисленными по начальному приближению управляющих параметров и полученными путем решения обратной задачи для одного из тестовых примеров. В данном случае размерность равномерной сетки составляла 21x21x5. Модель соответствовала обращенному 5-точечному элементу заводнения с 4 добывающими скважинами, расположенными в углах области фильтрации, и одной нагнетательной скважиной в центре. Добывающие скважины работали при постоянном дебите жидкости, а нагнетательная скважина производила полное восполнение объема отбора закачкой. В качестве управляющих параметров использовались 2 радиуса вариограмм в горизонтальной плоскости и 1 - в вертикальной, а также коэффициенты корреляционной зависимости. На рис. 9 и 10 вычисленные динамики значительно отличаются от «измеренных» для начального приближения управляющих параметров, но в точности совпадают с ними при использовании результатов решения обратной задачи.
а
« 182
* ♦
* * ♦ ★
-к * * замеры
начальное приближение адаптация
Рис. 9. Динамики забойного давления для одной из добывающих скважин в модели (начальное приближение соответствует завышенным значениям диапазонов вариограмм)
Рис. 10. Динамики забойного давления для
нагнетательной скважины (начальное приближение соответствует завышенным значениям диапазонов вариограмм)
Резюмируя, можно утверждать, что на сегодняшний день создан программный продукт, способный в автоматизированном режиме решать обратные задачи по уточнению коллекторских и иных свойств пласта. Это традиционное определение
множителей пористости и проницаемости по трем координатным направлениям, решение двух описанных задач геометрической идентификации, определение параметров геостатистической модели. Возможно уточнение кривых ОФП по регионам, определение параметров водоносного горизонта модели Фетковича. В разработке находится модуль определения динамики скин-фактора скважины в процессе ее эксплуатации. Количество достоверно определяемых параметров велико и может быть расширено по желанию пользователя пакета прикладных задач.
Статья написана в рамках выполнения государственного задания (тема «Научное обоснование новых экологически чистых технологий разработки месторождений углеводородов в сложных горно-геологических условиях на основе 3D-компьютерных экспериментов», № АААА-А16-116022510270-1).
ЛИТЕРАТУРА
1. Закиров С.Н., Коротаев Ю.П., Кондрат Р.М., Турниер В.Н., Шмыгля О.П. Теория водонапорного режима месторождений природного газа. М.: Недра, 1976. 240 с.
2. Coats K.H., Dempsey J.R., Henderson J.H. A New Technique Determining Reservoir Description from Field Performance Data // SPEJ. 1970. № 1. 9 p.
3. Абасов М.Т., Закиров И.С., Палатник БМ.Идентификация функций относительных фазовых проницаемостей при двухфазной фильтрации // ДАН СССР. 1990. Т. 312, № 4. С. 930-933.
4. Абасов М.Т., Закиров С.Н., Палатник Б.М. Адаптация геолого-математической модели газовой залежи при водонапорном режиме // ДАН СССР. 1989. Т. 38, № 2. С. 321-324.
5. Закиров С.Н., Васильев В.И., Гутников А.И., Коршунова Л.Г., Колбиков С.В. Прогнозирование и регулирование разработки газовых месторождений. М.: Недра, 1987. 295 с.
6. Закиров С.Н., Колбиков С.В., Палатник Б.М. Комплексные адаптирующиеся геолого-промысловые модели разработки газовых месторождений // Труды МИНГ им. И.М. Губкина. М., 1989. Вып. 214. С. 85-98.
7. Марчук Г.И., Агошков В.И., Шутяев В.П. Сопряженные уравнения и методы возмущений в нелинейных задачах математической физики. М.: Наука, 1993. 223 с.
8. Palatnik B., Aanonsen S.I., Zakirov I., Zakirov E. New Technique to Improve the Efficiency of History Matching of Full-Field Models // Proceedings of 4th European Conference on the Mathematics of Oil Recovery. Roros, Norway, 7-10 Ju^ 1994. 10 p.
9. Palatnik B.M., Zakirov I.S., Haugen S.A., van Roosmalen J.J. New Approach to Multiphase History Matching // Proceedings of the 7th European IOR Conference. Moscow, Russia, 27-29 October 1993. 10 p.
10. Закиров И.С. Развитие теории и практики разработки нефтяных месторождений. М.-Ижевск: Институт компьютерных исследований, НИЦ «Регулярная и хаотическая динамика». 2006. 356 c.
11. Закиров Э.С. Трехмерные многофазные задачи прогнозирования, анализа и регулирования разработки месторождений нефти и газа. М.: Изд. дом «Грааль», 2001. 302 с.
12. Nghiem L., Roson B. A Unified and Flexible Approach for Handling and Solving Large Systems of Equations in Reservoir Simulation // Proceedings of First and Second Forum on Reservoir Simulation, Alpbach, Austria, 1989. P. 501-550.
13. Азиз Х., Сеттари Э. Математическое моделирование пластовых систем. М.: Недра, 1982. 407 с.
14. Гилл Ф., Мюррей У., РайтМ. Практическая оптимизация. М.: Мир, 1985.509 с.
15. Yang P.H., Watson A.T. Automatic History Matching With Variable-Metric Methods // Paper SPE 16977 was first prepared for presentation at the SPE Annual Technical Conference and Exhibition. Dallas, 27-30 September 1987. SPERE, 1988. August. P. 995-1001.
16. Закиров Э.С. Upscaling в 3D компьютерном моделировании. М.: Изд-во «Книга и Бизнес», 2007. 344 c.
17. Закиров С.Н., Закиров Э.С., Индрупский И.М., Левченко В.С., Брадулина О.В., Цаган-Манджиев Т.Н. Вертикальное и 3D гидропрослушивание продуктивных пластов // Новые технологии освоения и разработки трудноизвлекаемых запасов нефти и газа и повышения нефтегазоотдачи: Сб. тр. VII Междунар. технологич. симпозиума. М., 2008. С. 49-63.
18. Zakirov E.S., Zakirov I.S. Aquifer Configuration Estimation Through Inverse Problem Solution // Paper SPE 51926 prepared for presentation at the 15th Reservoir Simulation Symposium. Houston, Texas. 14-17 February 1999. 2 р.
19. Закиров Э.С. Идентификация размеров и конфигурации водоносного пласта по данным разработки залежи нефти или газа // Наука и технология углеводородов. 2001. № 2. С. 71.
20. Zakirov I.S., Astl A., Zakirov E.S., Schweng K. History Matching for Lauchstaedt Underground Gas Storage // Paper SPE 39994 prepared for presentation at the Gas Technology Symposium. Calgary, Alberta, Canada. 15-18 March 1998. 10 p.
21. Закиров И.С., Hauenherm W., Закиров Э.С., Zipper H. History matching подземного хранилища Lauchstaedt // Газовая промышленность. 1997. № 10. С. 50-53.
22. Zakirov E.S., Indrupskiy I.M., Lubimova O.V., Shiriaev I.M. Geostatistically-consistent History Matching. // Proceedings of ECMOR XIV - 14th European Conference on the Mathematics of Oil Recovery. Catania, Italy. 8-11 September 2014. Vol. 2. P. 902-914.
23. Zakirov E.S., Indrupskiy I.M., Shiryaev I.M., Lyubimova O.V., Anikeev D.P. Advanced Geologically-consistent History Matching and Uncertainty Evaluation // Proceedings of ECMOR XV - 15th European Conference on the Mathematics of Oil Recovery. Amsterdam, the Netherlands. 29 August-1 September 2016. 26 p.
24. Закиров Э.С., Индрупский И.М., Любимова О.В., Ширяев И.М., Аникеев Д.П. Согласованная адаптация геостатистических моделей залежей нефти и газа // Доклады Академии наук. 2017. Т. 476. № 4. С. 421-425.
25. Goovaerts P. Geostatistics for Natural Resources Evaluation. New York: Oxford University Press, 1997. 483 р.
