Optimal number of Monte-Carlo simulations for hydrocarbon resources probabilistic estimation

Full Text

Введение

Симуляция Монте-Карло – широко используемая группа методов, в которых процесс или функция изучается путём многократных экспериментов (реализаций) с использованием определенной математической модели [1]. Симуляция используется в случаях, когда переменные функции неизвестны, но известно их распределение, а процесс или функция являются достаточно сложными для поиска аналитического решения.

Например, произведение двух случайных величин X, Y даёт третью случайную величину X·Y, дисперсия которой описывается достаточно громоздким уравнением. При использовании усечённых распределений определение дисперсии результата значительно усложняется.

Var(X·Y)=Var(X)·Var(Y) +

+ E(X)²·Var(Y) + E(Y)²·Var(X) (1)

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

Симуляция методом Монте-Карло позволяет рассчитывать процесс с достаточно большим количеством переменных. Основных недостатком метода является приближённость решения. Целью данной работы является определение оптимального числа симуляций при вероятностном подсчёте запасов для получения решения с необходимой погрешностью.

Методология

Симуляция Монте-Карло

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

y = f(X₁, X₂,…,X_m) (2)

Каждая переменная может быть представлена в виде случайной величины X_i, которая описывается плотностной функцией вероятности PDF_i (X_i) и вероятностной функцией распределения CDF_i (X_i) [2]. При этом значение случайной величины при определённой вероятности p рассчитывается обратной функцией от функции распределения:

X_i = CDF_i^-1(p) (3)

При расчёте реализации для каждой случайной величины генерируется случайное число, в результате получается m случайных чисел rand. В таком случае уравнение (2) приобретает вид (4):

$y=f\left(\begin{array}{c}CD{F}_1^{-1}\left(ran{d}_1\right), CD{F}_2^{-1}\left(ran{d}_2\right),... ,\\ CD{F}_m^{-1}\left(ran{d}_m\right)\end{array}\right)$ (4)

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

Таким образом, по результатам симуляции Монте-Карло получают массив результатов {y₁,y₂,…,y_n}, являющихся n-выборкой из функции, описывающей процесс. Предположив, что реализации между собой равновероятностны, итоговую функцию можно приближенно описать следующими статистическими величинами (табл. 1).

Таблица 1. Статистические величины для описания функции распределения

Table 1. Statistic values describing distribution function

Латинский гиперкуб и случайная выборка

Случайная выборка (Random sampler) – это способ получения случайной величины, при котором для его расчёта используется случайное число, определенное независимо от предыдущих реализаций. Однако, учитывая, что генератор случайных чисел в большинстве случаев не даёт равномерного распределения, случайная выборка может привести к значительному искажению итоговых результатов. Для преодоления влияния генератора случайных чисел используются разнообразные способы выборки, одним из которых является латинский гиперкуб.

Латинский гиперкуб (Latin-hypercube sampler) – это способ выбора значений из нескольких распределений, который обеспечивает более полный охват диапазона входных данных [3]. Это достигается разделением вертикальной оси накопленной функции распределения на n равных, непересекающихся интервалов (рис. 1). Таким образом, на горизонтальной оси также получается n интервалов с одинаковой вероятностью наступления, но разной шириной.

Следующим этапом является выбор значений внутри n интервалов. Для этого внутри каждого вертикального интервала случайной выборкой подбирается вероятность p_i (5), которой соответствует некоторое значение внутри интервала на горизонтальной оси (3).

p_i = p_{i min} + rand(p_{i max} – p_{i min}) (5)

где p_{i min}, p_{i max} – минимальное и максимальное значения вероятности i-го интервала.

Рисунок 1. Разбиение накопленной функции распределения на n интервалов

Figure 1. Partition of cumulative distribution function into n intervals

Виды распределений

Равномерное распределение – это вид распределения, при котором случайная величина равновероятностно может принять любое значение из диапазона от минимального до максимального значения. Для непрерывной случайной величины функция плотности равномерного распределения есть постоянная. Непрерывное равномерное распределение задается двумя параметрами: минимальное min и максимальное max значения случайной величины [4] (табл. 2).

Таблица 2. Основные виды распределений

Table 2. Main types of distributions

В случае вероятностного подсчёта ресурсов равномерное распределение используется для определения положения контакта в тех случаях, когда отсутствует какая-либо информация о положении контакта.

Нормальное распределение – наиболее часто используемый вид распределения, в котором функция плотности вероятности задаётся функцией Гаусса (табл. 2). Одной из отличительных особенностей данного распределения является его симметричность относительно среднего значения [5].

При вероятностном подсчёте ресурсов среднюю пористость и насыщенность пород часто описывают с использованием нормального распределения. Кроме того, в процессе вероятностного подсчёта ресурсов используют логнормальное распределение и гамма-распределение.

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

Усечение распределения – это преобразование распределения таким образом, что случайная величина не может выпадать за указанные рамки, т.е. случайная величина ограничена минимальным и максимальным значением (6).

$CD{F}_{tr}\left(x\right)=\frac{CDF\left(x\right)-CDF\left(min\right)}{CDF\left(max\right)-CDF\left(min\right)}$ (6)

Вероятностный подсчёт ресурсов

Структура X

Структура X представлена брахиантиклинальной складкой и была выявлена по данным 2D сейсмики. Нефтегазоносность региона подтверждена наличием месторождений поблизости. По данным месторождений-аналогов, а также учитывая глубины, предполагается, что в случае обнаружения залежи пластовым флюидом будет нефть. Подсчётные параметры были определены на основе данных месторождений-аналогов.

Учитывая высокую степень неопределённости, подсчёт ресурсов ведется вероятностным способом. Исходные данные представлены в табл. 3.

Таблица 3. Исходные данные для вероятностного подсчёта ресурсов

Table 3. Input data for the probabilistic resources estimation

Объёмный метод подсчёта ресурсов / запасов

При оценке ресурсного потенциала структуры был использован объёмный метод подсчёта ресурсов нефти. При расчёте геологических m_geo и извлекаемых m_rec ресурсов нефти производится многократное произведение общего объёма залежи BV на различные коэффициенты (7–8). Геологические VG_geo и извлекаемые VG_rec ресурсы попутного газа также рассчитываются путём перемножения соответствующих ресурсов нефти на газосодержание (9–10). Каждый из коэффициентов является независимой переменной, которую можно представить в виде непрерывной функции распределения.

$m_{geo}=BV·NTG·\phi ·S_o·FVF·{\rho }_o$ (7)

$m_{rec}=m_{geo}·RF$ (8)

$V{G}_{geo}=m_{geo}·R_S$ (9)

$V{G}_{rec}=m_{rec}·R_S$ (10)

Общий объём залежи. Различают два основных типа залежей: массивные и пластовые. Массивной называется та залежь, в которой непроницаемой является только покрышка. Внутренние пропластки внутри такой залежи гидродинамически связаны между собой.

Общий объём залежи BV зависит от геометрии пласта, его толщины H и положения контакта углеводорода с водой WC. Кровля пласта-коллектора структуры X была определена по данным геофизических исследований. Для простоты расчёта неопределённость, связанная с точностью картирования кровли пласта, не учитывается.

Учитывая степень заполнения ловушки месторождений-аналогов, предполагается наиболее вероятным заполнение залежи на половину. Учитывая замыкающую изолинию пласта (spill point), положение контакта описывается усечённым нормальным распределением (11). Толщина пласта была определена по данным гидродинамических исследований на близлежащих месторождениях и также описывается усечённым нормальным распределением (12).

$$\begin{gathered} WC=Trunc. Norm. Inv \\ ({\mu }_{WC}, {\delta }_{WC}^2, W{C}_{min}, W{C}_{max}) \end{gathered}$$ (10)

$$\begin{gathered} H=Trunc. Norm. Inv. \\ ({\mu }_H, {\delta }_H^2, H_{min}, H_{max}) \end{gathered}$$ (11)

Было проведено 100000 симуляций определения общего объёма залежи, в результате была получена гистограмма, представленная на рис. 2. Среднее значение полученной величины равно μ_BV = 635,6 млн м³, стандартное отклонение σ_BV = 400,0 млн м³.

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

Рисунок 2. Гистограмма распределения общего объёма залежи BV

Figure 2. Distribution histogram of the BV deposit total volume

Рисунок 3. Аппроксимация общего объёма залежи BV различными функциями распределений

Figure 3. Approximation of BV total bulk volume distribution by various distribution functions

Аналитический метод вероятностного подсчёта ресурсов

При подсчёте ресурсов объёмным методом производится перемножение переменных между собой. Для аналитического расчёта используется произведение двух независимых случайных величин. Пусть имеются две случайные величины X и Y со средним значением μ_X и μ_Y и дисперсиями Var(X) и Var(Y). Тогда произведение этих величин будет являться также случайной величиной со средним значением и дисперсией показанными в уравнениях (1) и (12).

${\mu }_{X·Y}={\mu }_X·{\mu }_Y$ (12)

Далее последовательно перемножается каждая следующая переменная на результат перемножения предыдущих переменных. Вначале рассчитывается эффективный объём залежи (13, 14). После рассчитывается поровый (15, 16) и нефтенасыщенный объёмы залежи (17, 18). Полученный нефтенасыщенный объём залежи перемножается на коэффициент усадки нефти, что даёт геологические ресурсы нефти в объёмных единицах (19, 20). Перемножив величину на плотность нефти, получаются геологические ресурсы нефти в массовых единицах (21, 22). Извлекаемые ресурсы равны произведению геологических ресурсов нефти и коэффициента извлечения нефти (23, 24). И, наконец, геологические и извлекаемые ресурсы попутного газа равны произведению соответствующих ресурсов нефти на газосодержание (25, 26).

${\mu }_{NV}={\mu }_{BV}·{\mu }_{NTG}$ (13)

${\delta }_{NV}^2={\delta }_{BV}^2·{\delta }_{NTG}^2+{\mu }_{BV}^2·{\delta }_{NTG}^2+{\mu }_{NTG}^2·{\delta }_{BV}^2$ (14)

${\mu }_{PV}={\mu }_{NV}·{\mu }_{por}$ (15)

${\delta }_{PV}^2={\delta }_{NV}^2·{\delta }_{por}^2+{\mu }_{NV}^2·{\delta }_{por}^2+{\mu }_{por}^2·{\delta }_{NV}^2$ (16)

${\mu }_{HCPV}={\mu }_{PV}·{\mu }_{So}$ (17)

${\delta }_{HCPV}^2={\delta }_{PV}^2·{\delta }_{So}^2+{\mu }_{PV}^2·{\delta }_{So}^2+{\mu }_{So}^2·{\delta }_{PV}^2$ (18)

${\mu }_{Qgeo}={\mu }_{HCPV}·{\mu }_{FVF}$ (19)

${\delta }_{Qgeo}^2={\delta }_{HCPV}^2·{\delta }_{FVF}^2+{\mu }_{HCPV}^2·{\delta }_{FVF}^2+{\mu }_{FVF}^2·{\delta }_{HCPV}^2$ (20)

${\mu }_{mgeo}={\mu }_{Qgeo}·{\mu }_{\rho o}$ (21)

${\delta }_{mgeo}^2={\delta }_{Qgeo}^2·{\delta }_{\rho o}^2+{\mu }_{Qgeo}^2·{\delta }_{\rho o}^2+{\mu }_{\rho o}^2·{\delta }_{Qgeo}^2$ (22)

${\mu }_{VGgeo}={\mu }_{mgeo}·{\mu }_{Rs}$ (23)

${\delta }_{VGgeo}^2={\delta }_{mgeo}^2·{\delta }_{Rs}^2+{\mu }_{mgeo}^2·{\delta }_{Rs}^2+{\mu }_{Rs}^2·{\delta }_{mgeo}^2$ (24)

${\mu }_{VGrec}={\mu }_{mrec}·{\mu }_{Rs}$ (25)

${\delta }_{VGrec}^2={\delta }_{mrec}^2·{\delta }_{Rs}^2+{\mu }_{mrec}^2·{\delta }_{Rs}^2+{\mu }_{Rs}^2·{\delta }_{mrec}^2$ (26)

Метод Монте-Карло вероятностного подсчёта ресурсов

Для симуляции вероятностного подсчёта ресурсов методом Монте-Карло генерируется двумерный массив случайных чисел размером {m, n}, где m – число переменных, а n – число симуляций. Затем, использовав случайное число массива в качестве значения накопленной функции вероятностного распределения CDF, получается уравнение (27), решив которое, можно получить переменную. Таким образом, решается обратная задача – нахождение аргумента, зная значение функции.

Например, параметр общего объёма залежи BV зависит от положения контакта и толщины пласта, которые в свою очередь являются переменными, моделируемыми с помощью нормального распределения. Квантиль, или значение переменной с накопленной вероятностью p, может быть выражен через обратную функцию ошибки erf ^-1 (p). Обратная функции ошибки может быть рассчитана путём разложения в ряд.

$F_X\left(x\right)={\int }_{-\infty }^x{f}_x\left(t\right)dt=Rand(0,1)$ (27)

где f_X – функция плотности распределения случайной величины X.

$X\left(p\right)=\mu +\delta \sqrt{2}·er{f}^{-1}(2p-1)$ (28)

Нахождение квантилей (значения переменной) гамма-распределения производится методом подбора: принимаются две точки l и r, значение накопленной фукнции в которых для l меньше p, а для r больше p (CDF(k,θ,l) < тp; CDF(k,θ,r) < p). Затем путём итеративного приближения друг к другу находится такое значение x, для которого выполняется условие |CDF(k,θ,x) – – p| ≤ eps, где eps – заданная точность переменной.

Проведя анализ над полученной выборкой, используя формулы из табл. 1, можно получить данные о статистическом распределении функции. Перцентили P10, P50, P90 рассчитываются либо по уравнению (3), используя вероятности 0,9, 0,5 и 0,1, либо путём ранжирования и разбиения множества полученных результатов.

Результаты

Случайная выборка или латинский гиперкуб

Для сравнения методов выборки случайных чисел была выбрана простая функция произведения двух случайных величин вида Y = X₁·X₂. Было проведено 20 расчётов с 10, 100, 1000 и 10000 симуляций, с использованием случайной выборки и с помощью латинского гиперкуба. Затем были рассчитаны ошибки среднего и среднеквадратичного отклонения относительно результатов, рассчитанных аналитическим способом. Результаты расчёта представлены на рис. 4. Как видно из рисунка, латинский гиперкуб во всех вариантах даёт меньшую погрешность, по сравнению со случайной выборкой (около 2 раз). Однако необходимо учитывать, что данный результат связан с несовершенством генератора случайных чисел. В настоящее время большинство генераторов случайных чисел завязаны на алгоритме псевдослучайных чисел, которые, в свою очередь, зависят от текущего такта центрального процессора компьютера.

Аналитическое решение

Формулы для аналитического решения представлены в уравнениях (13–26).

Однако данные уравнения были модифицированы, т.к. в качестве входных данных использовались усечённые распределения. В табл. 4 представлены результаты аналитических расчётов. Полученные результаты использовались для сравнения с результатами, полученными путём симуляции Монте-Карло при разном количестве симуляций.

Таблица 4. Результаты аналитического вероятностного подсчёта ресурсов

Table 4. Results of analytical probabilistic resources estimation

Число симуляций

Были проведены серии расчётов аналитическим методом и методом имитационного моделирования для различного числа симуляций. Для исключения влияния генератора случайных чисел расчёт повторялся 50 раз, каждый раз генерировались новые массивы случайных чисел. Анализировалось среднее значение погрешности при 10, 50, 100, 200, 500, 1000, 10000, 50000 и 100000 симуляций подсчёта ресурсов (табл. 5).

Оценивалась погрешность результатов относительно аналитического решения, которое даёт точный ответ. Увеличение числа симуляций повышает точность вычислений. Кроме того, расчёты чувствительны к генератору случайных чисел, из-за чего возможны случаи, когда большое число симуляций может снизить точность. Эмпирическим путём были получены формулы, с помощью которых можно определить требуемое число симуляций (рис. 5).

Таблица 5. Влияние количества симуляций на ошибку при симуляции Монте-Карло

Table 5. Effect of the number of simulations on the Monte Carlo simulation error

Рисунок 4. Сравнение методов выборки случайного числа

Figure 4. Comparison of random number sampling methods

Рисунок 5. Средняя погрешность при разном количестве симуляций

Figure 5. Average error for different number of simulations

Выводы

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

Симуляция Монте-Карло позволяет приближенно находить решение для сложных функций, в которых может быть большое количество переменных. Точность решения зависит от количества симуляций, однако после определенного порога дополнительное увеличение количества симуляций практически не увеличивает точность расчёта.

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

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

ДОПОЛНИТЕЛЬНО

Источник финансирования. Автор заявляет об отсутствии внешнего финансирования при проведении исследования.

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

ADDITIONAL INFORMATION

Funding source. This study was not supported by any external sources of funding.

Competing interests. Author declare that he has no competing interests.

About the authors

Raman M. Sadykov

Author for correspondence.
Email: r.sadykov@niikmg.kz
ORCID iD: 0000-0002-5936-2036
Kazakhstan, Astana

References

Supplementary files

Supplementary Files

Action

1. JATS XML

Download

2. Figure 1. Partition of cumulative distribution function into n intervals

Download (118KB)

Indexing metadata

3. Figure 2. Distribution histogram of the BV deposit total volume

Download (140KB)

Indexing metadata

4. Figure 3. Approximation of BV total bulk volume distribution by various distribution functions

Download (137KB)

Indexing metadata

5. Figure 4. Comparison of random number sampling methods

Download (93KB)

Indexing metadata

6. Figure 5. Average error for different number of simulations

Download (92KB)

Indexing metadata