CC BY
76
Математическое моделирование
УДК 519.218.7, 519.218.8
Спектральные методы стохастического моделирования гауссовских процессов в задачах автоадаптации
И. Р. Минниахметов, Д. А. Митрушкин
Отдел 11, сектор 3 Институт прикладной математики им. М.В. Келдыша РАН Миусская пл., д.4, Москва, 125047, Россия
В данной работе рассматривается задача адаптации гидродинамических моделей углеводородных месторождений на основе сопоставления гидродинамических расчётов с историческими данными наблюдений. В качестве варьируемых физических свойств пласта используются фильтрационно-ёмкостные свойства (ФЕС) среды: пористость и проницаемость. Ключевым моментом процесса адаптации является выбор способа параметризации полей ФЕС. В работе предложен эффективный метод спектральной параметризации полей ФЕС на основе алгоритма разложения Холецкого матрицы ковариации условного процесса в Фурье-пространстве. Данный подход позволяет существенно сократить число запусков гидродинамических расчётов. Проведён сравнительный анализ результатов расчётов предложенного алгоритма со стандартным спектральным методом на тестовой модели PUNQ-S3.
Ключевые слова: стационарные процессы, гауссовские процессы, спектральный метод, преобразование Фурье, разложение Холецкого, автоадаптация.
1. Введение
Математическое моделирование является неотъемлемой частью создания проектов разработки нефтегазовых месторождения. С помощью построенных геолого-гидродинамических моделей проводится расчёт прогнозных вариантов разработки, оценивается эффективность планируемых геолого-технологических мероприятий, направленных на увеличение добычи и коэффициента извлечения углеводородов, принимаются оперативные решения по эксплуатации месторождения. Качество гидродинамической модели оценивается по степени её адаптации к историческим данным разработки месторождения, а также соответствия параметров модели результатам экспериментов и исследований. Одним из инструментов адаптации является изменение полей фильтрационно-ёмкостных свойств (ФЕС) пласта: пористости и проницаемости. Однако ручная модификации ФЕС требует большого количества времени и высокой квалификации специалистов, поэтому актуальной является разработка автоматизированных алгоритмов подбора свойств модели.
В основе гидродинамической модели обычно лежит усреднённая геологическая модель, построенная на основе скважинных данных и результатов анализа сейсмических исследований и фациально-литологических свойств пород [1]. Полученные при этом значения ФЕС с достаточной степенью достоверности определены только вблизи скважин, а в межскважинном пространстве заданы весьма условно. Поэтому при проведении адаптации гидродинамической модели в качестве варьируемых параметров можно использовать значения полей ФЕС в межскважинной области, а сами поля рассматривать в виде реализаций стационарных случайных процессов [2]. При проведении адаптации важно не только обеспечить совпадение результатов моделирования с историей разработки, но и получить физически обоснованные поля ФЕС, обладающие свойством устойчивости по отношению к неточности задания исходных данных.
Статья поступила в редакцию 7 ноября 2012 г.
Работа выполнена при частичной финансовой поддержке Минобрнауки России (госконтракт № 07.514.11.4008) и РФФИ (проект № 12-01-00793-а).
Ввиду большого числа расчётных ячеек гидродинамической модели (Лл ~ 106) варьирование значений ФЕС во всех точках межскважинной области является трудноразрешимой задачей даже для вычислительной техники современного уровня. Отсюда возникает необходимость введения эффективной параметризации полей ФЕС в межскважинной области с числом параметров Хр ^ Л^. На сегодняшний день используются три основных типа параметризации [3]: региональная параметризация, метод пилотных точек, глобальная параметризация.
Региональная параметризация заключается в разбиении полей ФЕС на регионы «однородности», каждый из которых регулируется своим параметром. Отсутствие строго обоснования такого разбиения и получившиеся в результате оптимизации резкие переходы на границах регионов являются несомненными недостатками данного подхода. Тем не менее региональная параметризация является стандартным методом адаптации гидродинамических моделей.
Метод пилотных точек использует в качестве параметров значения полей пористости и проницаемости в конечном наборе точек, называемых пилотными. Значения ФЕС в межскважинном пространстве вычисляются на основе скважинных данных и значений в пилотных точках, с помощью геостатистических методов [2]: кригинг, последовательная гауссовская симуляция и др. Таким образом, данный подход позволяет провести адаптацию модели, сохранив геостатистические свойства полей ФЕС: корреляционную функцию, среднее значение и дисперсию. Недостатками метода пилотных точек является сложный алгоритм выбора количества и координат пилотных точек, а также локальные «выбросы» значений ФЕС в области пилотных точек.
В качестве параметров для глобального типа параметризации используются свойства модели, не зависящие от координат, например: средние значения и дисперсия полей ФЕС, радиусы и направления корреляции, форма проводящих каналов и др. Однако перечисленные параметры не позволяют локально варьировать значения ФЕС.
В данной работе предлагается использовать в качестве параметров амплитуды гармоник Фурье-образов полей ФЕС [4]. Использование спектрального подхода [5] позволяет не только провести адаптацию гидродинамической модели, но и значительно сократить количество запусков гидродинамического симулятора.
2. Постановка задачи
Для простоты изложения, не ограничивая общности, вместо функций пористости и проницаемости пласта будем рассматривать скалярную функцию филь-трационно-ёмкостного свойства X(г), где г — радиус вектор. Опишем входные данные процесса автоадаптации:
1. Средние значения ФЕС пласта а(г), полученные на основе геологических карт пористости и проницаемости.
2. Значения ФЕС (I) вдоль стволов скважин, где ] — индекс скважины, I — расстояние от устья скважины.
3. Данные истории разработки ¿^(Д^) на определённый период времени Д^, заданные с ошибкой о^к, где % — наблюдаемое значение, например, значение забойного давления, дебит газа, нефти и воды.
Геометрия пласта описывается трёхмерной сеткой П. Требуется построить поле ФЕС X(г) на сетке П, удовлетворяющее следующим условиям:
1. Поле X (г) совпадает со значениями отсчётов тт в ячейках, через которые проходят скважины:
X(г,т) = (щ (1))т ,т = 1,М, (1)
где — радиус вектор центра ячейки т сетки, через которую проходит скважина ], ()т — усреднение скважинных данных, попавших в ячейку т, М -количество ячеек сетки, через которые прошли скважины.
2. Сумма квадратов отклонений значений (Ак), рассчитанных на основе поля X(г), от наб. ния ошибки :
ля X(r), от наблюдаемых значений d°js(Atk) не превышает заданного значе-
^ щ
5& = ЕЕЕ [(^(Аь) - ^(X))/агок]2 < £, (2)
1=1 ]=1 к=1
где Ха — число наблюдаемых значений, Хг — число периодов разработки, Хт — число скважин.
Для уменьшения числа переменных будем использовать параметризацию поля X(r):
X (r) = А? (3)
где ? — параметры, А — оператор отображения.
В таком случае задача (2) преобразуется следующим образом:
S(r) = S(Ap). (4)
Однако данная задача является плохообсуловленной [3]. Во-первых, решение задачи не является единственным. Во-вторых, для конкретного решения задачи результаты расчётов симулятора могут практически не зависит от значений некоторых параметров. Таким образом данные параметры могут принимать практически любые значения.
В данной работе будем использовать баесовский подход, позволяющий уменьшить число обусловленности задачи. Согласно баесовскому формализму, данные наблюдения d°jS(Atк) являются случайной выборкой с параметрами ?, а выражение (4) представляет собой функцию правдоподобия f(dohs\р). Теперь положим, что априорное распределение д(р) существует, тогда параметры ?можно рассматривать как случайные величины с апостериорным распределением (АР):
p^cf(dobs\p)g(p) (5)
где — константа.
Предполагая, что распределения являются гауссовскими, выражение (5) примет вид:
?н cexp(- [S (?) + (? - ?Р)ТС^1(? - ?р)] /2) (6)
где jj,p — математическое ожидание вектора р, Ср — матрица ковариации вектора р.
Таким образом задача адаптации сводится к оценке апостериорного максимума (МАР) параметров р:
Рмар (dobs) = argmaxf(dobs\p)g(p), (7)
р
что эквивалентно
Рмар(dobs) = argmin [S(?) + (?- ?Р)ТС-1(?- ?р)] /(XdXtXw + Np), (8) р
и определению её неопределённости.
3. Параметризация
Решение поставленной задачи будем проводить в предположении [2], что ФЕС среды являются случайными процессами, представимыми в виде:
X (г) = а(г) + Ф(г) (9)
где Ф(г) — гауссовский стационарный случайный процесс [6] с нулевым средним и заданным корреляционным оператором В (г, г') = В(г — г'), зависящим только от разности векторов р = г — г'.
Будем считать, что корреляционная функции В(р) имеет гауссовский вид:
В(р) = а2 еМ—М2 /&)■ (10)
Как известно [6], гауссовский стационарный процесс может быть представлен в виде:
Ф(г) = ! д(й)<п{й)(13ш, (11)
где ш — гармоники Фурье-образа Ф(г), функция д(ш) такая что д(ш)д*(ш) = В(ш), В(ш) — Фурье-образ функции В(р), 'д(ш) — независимые стандартные нормальные случайные величины.
Сеточные значения процесса Ф(г) представляются в виде:
ф = ЕЩ (12)
где Е — матрица дискретного преобразования Фурье [7], Ю — диагональная матрица с элементами Бц = д(ш{), щ — гармоники Фурье-образа Ф.
Далее реализации процесса Ф будем обозначать 'ф. В работе [5] показано, что в случае, когда реализации -ф случайного процесса Ф удовлетворяют условию (1), компоненты вектора не являются независимыми и процесс моделирования реализаций ф значительно усложняется. Такие процессы будем называть условными и обозначать фсопЛ.
Стандартный подход [8] моделирования реализаций условного гауссовского процесса фсопЛ основан на коррекции безусловных гауссовских полей (КБГП) (12) посредством учёта невязок Д^ между значениями на скважинах (1) и значениями полученной реализации ф, где к — индекс значения ФЕС вдоль ствола скважин.
На первом этапе генерируется реализация ф, используя выражение (12). Далее определяются невязки Д^ и с помощью метода кригинг вычисляются значения невязок Д(г) во всех ячейка сетки П. Результирующая реализация представляется в виде:
ф)сопА = >ф + Д (13)
Согласно (12), матрица Ю отвечает за амплитуды гармоник Фурье-образа вектора 'фсопЛ. Так как Фурье-образ функции (10) является функцией Гаусса с шириной Кш ~ 2ъ/К, то наиболее значимые гармоники шзгдп находятся в шаре радиуса 3ВШ с центром в начале координат. Отсюда следует, что в качестве параметров для решения задачи (4) могут выступать компоненты вектора г] , соответствующие гармоникам ш8гдп. Однако, несмотря на то, что число значимых гармоник
много меньше количества ячеек сетки ^ Ил, количество параметров для
оптимизации остаётся достаточно высоким.
В данной работе мы предлагаем использовать подход [5], позволяющий свести процесс оптимизации к нескольким последовательным оптимизациям с малым числом параметров. Метод основан на алгоритме эффективного разложения
Холецкого (ЭРХ) [9] матрицы ковариации £ условного процесса фсопЛ в Фурье-пространстве. В этом случае реализации (13) представимы в виде:
ifcond = FD^ = FDLW £ (14)
где Lw — множитель ЭРХ матрицы ковариации £ в Фурье пространстве.
Другими словами, если для безусловного процесса Ф в качестве пространства с независимыми элементами достаточно использовать Фурье-пространство,
то для условных процессов фсопа требуется дополнительное отображение Lw. Заметим, что аналогично (12) в выражении (14) наиболее значимые компоненты £ определяются матрицей D.
Благодаря тому, что матрица Lw является нижней треугольной [10], значение амплитуды произвольной гармоники зависит только от значений компонент £i,£2,-£fc вектора £ :
A„k = f ), (15)
где АШк — амплитуда гармоники Wk, £,k — компоненты вектора £ .
Следовательно, процесс оптимизации задачи (4) может быть сведён к нескольким последовательным оптимизациям: сначала оптимизация по компонентам вектора £ , соответствующая низкочастотным гармоникам, далее по среднечастот-ным и высокочастотным, тем самым покрывая всю область гармоник us%gn. Таким образом размерность пространства параметров существенно уменьшается, что позволяет значительно сократить количество расчётов гидродинамического симулятора и ускорить процесс автоадаптации.
4. Оценка неопределённости
Помимо определения оптимальной реализации X(рмар) важнейшим этапом адаптации является оценка неопределённости. Для решения данной задачи необходимо построить АР (6). Однако единственным практическим способом получения АР является оценка параметров распределения по выборке. Известно [11], что используя метод Монте-Карло по схеме марковской цепи (МКМЦ) можно корректно построить необходимую выборку, однако данный подход требует большого количество запусков гидродинамического симулятора, что является трудо-затратной процедурой. В данной работе используется метод RML(Randomized Maximum Likelihood) [12]. Для линейных задач доказано [13], что выборка, полученная данным методом, входит в класс корректности метода МКМЦ. Однако для нелинейного случая строгое теоретическое обоснование отсутствует. Тем не менее на частных примерах даже на нелинейных задачах RML даёт [14] разумную характеристику неопределённости.
5. Тестовая модель
Продемонстрируем эффективность предложенного подхода на примере задачи PUNQ-S3 [3] , основанной на модели реального нефтегазового месторождения.
Месторождение размерами 3420 х 5039 х 25 метров имеет куполообразный вид с газовой шапкой, газонефтяным и водонефтяным контактами (рис. 1). На востоке и юге ограничено разломом, а со стороны севера и запада - водонапорным горизонтом. Разработка месторождения осуществляется 6 добывающими скважинами, расположенными по всему газонефтяному контакту. Период истории разработки составляет 16.5 лет. Для каждой скважины заданы временные интервалы их работы с соответствующими значениями забойного давления, газонефтяного фактора и величины обводнённости. Общая добыча нефти за весь период разработки составляет 3, 87 х 106м3.
Расчётная сетка содержит 19 х 28 х 5 ячеек, из которых 1761 являются активными. Исходная геолого-гидродинамическая модель содержит все необходимые свойства флюидов и водонапорного горизонта. Для каждой скважины заданы кривые пористости, проницаемости по вертикали и горизонтали.
Для реалистичности данных к скважинным значениям ФЕС добавлен гауссов-ский шум с дисперсией 15%. Данные по давлению и дебитам флюидов зашумлены с коррелированными во времени величинами, чтобы отразить систематический характер ошибок наблюдений. Уровень шума для давления закрытия выбирался в 3 раза меньше, чем для давления в пласте, соответственно, 1 бар и 3 бара, чтобы более точно отразить момент закрытия. Уровень шума для газового фактора установлен на уровне 10% до прорыва газа и 25% после прорыва газа, тем самым отражая разницу между растворенным и свободным состоянием газа. Аналогично для обводнённости 2% до и 5% после прорыва воды.
Рис. 1. Модель (вид сверху)
6. Результаты расчётов
В данном параграфе представлены результаты адаптации, используя подходы КБГП и РХ. Для начала, все компоненты вектора ] были отсортированы так, чтобы Фурье-образ корреляционной функции В (г) убывал. Таким образом компоненты вектора ] были упорядочены по степени их влияния на реализацию процесса (9). Далее компоненты вектора ] разбиты на Р групп, для каждой из которых проводилась оптимизация методом контролируемого случайного поиска [15] и ограниченным числом запуска симулятора (.
В таблице 1 представлены основные результаты оптимизации для различных значений Р, (.
Таблица 1
Результаты оптимизации и оценки неопределённости
Метод Количество итераций ( Число групп Р Ошибка адаптации 5 (р)
КБГП 500 10 1.35
ЭРХ 500 10 1.24
КБГП 1000 50 1.19
ЭРХ 1000 50 1.05
Также сравнивалось количество запусков симулятора при заданной точности расчёта. Для достижения точности е = 1 методу ЭРХ потребовалось 1083 запусков гидродинамического симулятора, в то время как методу КБГП — 1683. Такой выбор точности обусловлен тем фактом, что точность расчётов должна быть согласованна с точностью входных данных.
Наилучшие результаты адаптации для исследуемых методов представлены на рисунках 2, 3 и 4.
Рис. 2. Результаты адаптации для накопленной добычи нефти
Рис. 3. Результаты адаптации для накопленной добычи газа
Рис. 4. Результаты адаптации для накопленной добычи воды
Сопоставив полученные результаты с данными из работы [3] (рис.5), можно сделать вывод, что предложенный метод ЭРХ позволяет составить достаточно точный прогноз с приемлемой степенью неопределённости.
6
р-н
.....Истинное значение
TNO-1 TNO-2 TNO-a Amoco Amoco Elf IFP- IFP- NCC-. NCC- NCO C)PY iso aniso STM Oliver Oliver GA MCMC
iso
Рис. 5. Прогнозные значения накопленной добычи нефти за 16,5 лет разработки и неопределённость прогноза, полученные с использованием различных подходов автоадаптации
В работе предложены два спектральных подхода в применении к задаче адаптации месторождений. Использование метода ЭРХ позволяет не только совпасть с историей разработки, но и сохранить геостатистические свойства полей ФЕС. Методы адаптации опробованы на тестовой модели PUNQ-S3. Согласно полученным результатам, использование подхода РХ позволяет значительно сократить количество запусков гидродинамического симулятора, что является ключевым моментом для прогнозирования добычи углеводородов и оценки неопределённости полученных результатов.
1. Закревский К. Е. Геологическое 3D моделирование. — Москва: ООО ИПЦ «Маска», 2009. [Zakrevskiyj K. E. Geologicheskoe 3D modelirovanie. — Moskva: OOO IPC «Maska», 2009. ]
2. Dubrule O. Geostatistics in Petroleum Geology // Am. Assn. Petr. Geol. — 1998. — No 38.
3. Floris F. J. T. et al. Methods for Quantifying the Uncertainty of Production Forecasts; a Comparative Study // Petroleum Geoscience. — 2001. — Vol. 7.
4. Jafarpour B., McLaughlin D. B. Reservoir Characterization with the Discrete Cosine Transform // SPE J. — 2009. — Vol. 14, No 1. — Pp. 182-201.
5. Минниахметов И. Р. Стохастическое моделирование условных гауссов-ских процессов. — 2011. [Minniakhmetov I. R. Stokhasticheskoe modelirovanie uslovnihkh gaussovskikh processov. — 2011. ]
6. Крамер Г., Лидбеттер М. Стационарные случайные процессы. — М.: Мир, 1969. [Kramer G., Lidbetter M. Stacionarnihe sluchayjnihe processih. — M.: Mir,
7. Lidl R., Pilz G. Applied Abstract Algebra. — New York: Wiley, 1999.
8. Dietrich C. R. Computationally Efficient Generation of Gaussian Conditional Simulations over Regular Sample Grids // Math. Geol. — 1993. — Vol. 25, No 4. — Pp. 439-451.
9. Haugh M. The Monte Carlo Framework, Examples from Finance and Generating Correlated Random Variables // IEOR E4703. — 2004.
10. Bau III D., Trefethen L. Numerical linear algebra. — Philadelphia: Society for Industrial and Applied Mathematics, 1997. — Pp. 172-180.
11. Oliver D., Cunha L., Reynolds A. Markov Chain Monte Carlo Methods for Conditioning a Permeability Field to Pressure Data // Math. Geol. — 1997. — Vol. 29.
7. Заключение
Литература
1969.]
12. Feng T., Mannseth T., Aanonsen S. I. Randomized Maximum Likelihood with Permeability Samples Generated by a Predictor Corrector Technique // University of Bergen 2009. Society of Petroleum Engineers. — 2012.
13. Kitanidis P. Quasi-Linear Geostatistical Theory for Inversing // Water Resour. Res. — 1995. — Vol. 31. — P. 2411-2419.
14. Liu N., Oliver D. S. Evaluation of Monte Carlo Methods for Assessing Uncertainty // SPEJ. — 2003. — Vol. 8, No 2. — Pp. 1-15.
15. Kaelo P., Ali M. M. Some Variants of the Controlled Random Search Algorithm for Global Optimization // J. Optim. Theory Appl. — 2006. — Vol. 130, No 2. — Pp. 253-264.
UDC 519.218.7, 519.218.8 Spectral Approaches of Gaussian Conditional Simulations in History Matching Problems
I. R. Minniakhmetov, D. A. Mitrushkin
Departament 3, group 11 Moscow Institute of Physics and Technology 1a, Kerchenskaya str., Moscow, 117303, Russia
In the paper we consider optimization problem arised in history matching of reservoir model. In such type of problems the unknown parameter often is a distributed field of the physical quantity such as permeability and porosity fields. The way the parameter field is parametrized greatly influences efficiency of the complete optimization approach. We propose an efficient technique of a spectral-domain parameterization based on the Cholesky decomposition of the covariance matrix in Fourier domain. The approach significantly reduce the number of simulation runs. A comparative analysis of history matching of the proposed algorithm and standard spectral method is performed using PUNQ-S3 test model.
Key words and phrases: stationary process, Gaussian fields, spectral approach, Fourier transformation, Cholesky decomposition, history matching.
