Анализ кривой падения добычи нефтяных и газовых скважин

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

Ключевые моменты практического анализа падения добычи

Традиционный подход к прогнозированию добычи, имеющий почти столетнюю историю, — это анализ кривой падения (DCA). Это эмпирический метод прогнозирования производительности нефтяных и газовых скважин. Неопределенности этого подхода зависят от сложности будущего поведения скважины и умения аналитика их прогнозировать. Методики DCA позволяют построить кривую производительности как функцию от времени. Некоторые резервуары, такие как нефтяные с сильным притоком воды, имеют достаточно энергии, чтобы поддерживать относительно постоянный расход нефти или газа в течение длительного периода времени. Как правило, производительность снижается со временем после завершения скважины и начала добычи. Целью DCA анализа является составление прогноза будущих темпов добычи и определение ожидаемых конечных извлекаемых запасов.

Производительность имеет начальный переходный период потока (initial transient flow period), за которым следует период потока с преобладанием границ (boundary-dominated flow period). В переходный период пластовое давление на границе потока остается постоянным при начальном пластовом давлении, а граница потока перемещается наружу из скважины через пласт. Эта часть скважины характеризуется очень высокими темпами падения. Когда граница потока достигает фактической границы пласта или встречается с границей потока другой скважины, давление пласта начинает падать, и скважина входит в период течения с преобладанием границы. Именно в этот период можно использовать традиционные методы падения (т. е. методы Arps).

С развитием методов машинного обучения (ML), применяемых в нефтегазовой отрасли некоторые алгоритмы превосходят традиционные методы DCA, такие как Arps, Duong и stretched exponential production decline (SEPD). Аналогичное улучшение наблюдалось в применении алгоритмов ML для типовых прогнозов скважин.

Arps модель

Теория анализа всех кривых спада начинается с концепции номинальной (мгновенной) скорости падения (decline rate), которая определяется следующим уравнением

$a = - \frac{\Delta{q}/q}{\Delta t} = - \frac{1}{q} \frac{\Delta q}{\Delta t}$

Другой способ представления скорости снижения основан на расходе (rate) и константе показателя снижения (decline exponent constant) .

$a = k q^{b}$

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

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

Экспоненциальный спад (Exponential) — b = 0
Гиперболический спад (Hyperbolic) — b имеет значение, отличное от 0 или 1
Гармонический спад (Harmonic) — b = 1

В расчетах должны быть согласованы единицы измерения между расходом и скоростью снижения. Например, если расход указан в stb/day, скорость снижения должна быть в единицах 1/day, а время должно быть в единицах day. В этом случае, если время указано в месяцах, нужно будет добавить количество дней, соответствующее конкретному рассматриваемому месяцу (28, 29, 30 или 31 день) на каждой итерации.

Для приведенного выше уравнения есть два особых случая: когда b = 0 и когда b = 1 .
Случай b=0 известен как экспоненциальный спад и вызовет ошибку деления на ноль в приведенном выше уравнении. Применяя теорию пределов, когда стремится к 0, уравнение становится:

$q(t)=q_i e^{-a_i t}$

Случай b=1 известен как гармонический спад. Он выражается как особый случай, поскольку исходная теория не предсказывает коэффициенты больше 1. Однако уравнение Arps можно применять для переходных потоков (transient flows), типичных для нетрадиционных скважин. Таким образом, коэффициенты , равные 2,5 или более, являются обычными. В любом случае, уравнение гармонического спада имеет вид:

$q(t)=\frac{q_i}{1+a_i t}$

Более подробный вывод формул Arps метода можно посмотреть в статье Decline Analysis Theory

Альтернативные модели DCA

Более подробно типы DCA анализов можно ознакомиться в PetroWiki
Модель Stretched Exponential Production Decline (SEPD) предлагается для того, чтобы адекватно моделировать показатели падения добычи из нетрадиционной залежи. SEPD основана на идее, что несколько спадающих систем составляют одну. Скорость добычи снижается со временем в соответствии со следующим соотношением:

$q(t)=q_iexp[-(\frac{t}{\tau})^n]$

Где:
q_i — начальная скорость
$\tau$ — характерный параметр времени
— параметр экспоненты
— период времени (обычно дни или месяцы).

Метод Duong был разработан специально для нетрадиционных коллекторов с очень низкой проницаемостью. Форма кривой подходит для скважин, которые демонстрируют длительные периоды неустановившегося потока. Метод Duong достигает конечного значения eur (expected ultimate recoverable reserves — ожидаемые конечные извлекаемые запасы). В этом случае предполагается увеличение плотности трещин с течением времени из-за истощения давления.
Уравнение падения Duong

$q=q_1t^{-m}exp[\frac{a}{1-m}(t^{1-m}-1)]$

Где:
qi -начальная скорость падения,
a — это точка пересечения графика обратного материального баланса q/N_p в двойном логарифмическом масштабе в зависимости от времени,
— наклон графика обратного материального баланса q/N_p в двойном логарифмическом масштабе от времени, t
— период времени (обычно дни или месяцы)

Метод PowerLow основан на модели Arps. Единственное отличие состоит в том, что спад в модели Arps является постоянным, тогда как модель PowerLow рассматривает спад как функцию степенного закона.

$q(t)=q_iexp(-D_{\infty}t-D_it^n)$

Где
q_i — это точка пересечения скорости или q(t=0)
D_i — это константа снижения
$D_\infty$ — это константа снижения
— это временная экспонента
— это период времени (обычно дни или месяцы)

ARIMA (autoregressive integrated moving average) моделирует данные временных рядов для прогнозирования (т. е. для предсказания будущих точек ряда) таким образом, что происходит следующее:

p — учитывается закономерность роста/снижения данных (авторегрессионный компонент);
d -учитывается скорость изменения роста/снижения данных (интегрированный компонент);
q -учитывается шум между последовательными временными точками (компонент скользящего среднего). Например, p=2, d=1 и q=0

from statsmodels.tsa.arima.model import ARIMA  model_AR = ARIMA(ts_log, order=(2, 1, 0)) results_ARIMA_AR = model_AR.fit() plt.figure(figsize=(10, 5)) plt.plot(ts_log_diff_active, color='blue', label='Original') plt.plot(results_ARIMA_AR.fittedvalues, color="red", label='Fitted') plt.title("RSS: %.3f" % sum((results_ARIMA_AR.fittedvalues - ts_log_diff_active) ** 2)) plt.legend() plt.show()

Подробнее с ARIMA можно познакомиться в работе «Machine Learning in the
Oil and Gas Industry»

Реализация Arps метода на Python

import numpy as np import matplotlib.pyplot as plt  def arps(t, qi, ai, b = 0):     t = np.array(t)     if b > 0:         # hyperbolic cases         q = qi*np.power(1 + b*ai*t, -1.0/b)     elif b == 0:         # exponential case         q = qi*np.exp(-ai*t)         # harmonic cases     elif b == 1:         q = qi / (1 + ai*t)     else:         # invalid b factor         raise Exception('Invlid b factor')     return q  qi = 1000 ai = 0.01 t = [x for x in range(1, 3650)]  q_exp = arps(t, qi, ai, 0) q_hyp = arps(t, qi, ai, 0.5) q_har = arps(t, qi, ai, 1)  plt.plot(t, q_exp, '-', label='Exponential decline') plt.plot(t, q_hyp, '-', label='Hyperbolic decline') plt.plot(t, q_har, '-', label='Harmonic decline') plt.legend(loc='upper right') plt.title('DCA') plt.xlabel('Time') plt.ylabel('Production') plt.plot()

В данном примере взят период около 10 лет — 3650 дней. Начальная производительность скважины qi = 1000 stb/day. Начальная скорость падения выбрана ai = 0.01 1/day. Представлены три кривые спада: для случая b = 0 — экспоненциальный спад; b = 0.5 — гиперболический спад; b = 1 — гармонический спад.

Модифицированная гиперболическая модель (Modified Arps)

Использование гиперболических уравнений с экспоненциальными коэффициентами выше 1 ведет к переоценке срока службы скважины, поскольку дебит снижается очень медленно на позднем этапе ее эксплуатации. Модифицированная гиперболическая модель была создана для решения этой проблемы. Основная идея модифицированной гиперболической модели заключается в том, чтобы определить, когда снижение становится слишком малым, чтобы быть реалистичным, а затем переключить снижение с гиперболического на экспоненциальное. Существует три способа определения начальной скорости снижения для гиперболического снижения: номинальный спад (Nominal decline), эффективный спад по касательной (Tangent effective decline) и эффективный спад по секущей (Secant effective decline).

Экспоненциальный спад

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

$D = \frac{ln(\frac{q_i}{q})}{t}$

Эффективное снижение для определенного периода времени, обычно 1 год

$D_e = \frac{q_i - q}{q_i}$

Эффективное снижение как функция от номинального

$D_e = 1 - e^{-D}$

Номинальный спад как функция от эффективного

Гиперболический спад

Номинальный спад

$D_i = \frac{(\frac{q_i}{q})^b}{bt}$

Тангенциальный эффективный спад, где qi и q считываются с касательной линии

$D_{ei} = \frac{q_i - q}{q_i}$

Секущий эффективный спад, где qi и q считываются с секущей линии

$D_{esi} = \frac{q_i - q}{q_i}$

Номинальный спад как функция от тангенциального эффективного спада

$D_i = -ln(1 - D_{ei})$

Номинальный спад как функция от секущего эффективного спада

$D_i = \frac{(1 - D_{esi})^{-b} - 1}{b}, b \neq 0$

Где:
D — номинальная экспоненциальная скорость снижение, 1/time;
Di — начальная номинальная скорость снижения (t = 0), 1/time;
Dei — начальная эффективная скорость снижения с касательной линии, 1/time;
Desi — начальная эффективная скорость снижения с секущей линии, 1/time;
qi — мгновенная производительность в начальный момент времени;
q — мгновенная производительность в момент времени t;
t — время;
e — основание натурального логарифма (2.718…);
b — гиперболический показатель (описывает, как начальная скорость снижения, Di, изменяется со временем; обычно варьируется от 0 до 1).

Реализация Modified Arps на Python

import numpy as np import matplotlib.pyplot as plt  def modified_arps(t, qi, di, b, lim = 0.05, type = 'Secant'):     if type == 'Secant':         # Nominal decline as a function of secant effective decline         d_nom = (np.power(1 - lim, -b) - 1.0) / 365/ b     elif type == 'Tangent':         # Nominal decline as a function of tangent effective decline         d_nom = -np.log(1 - lim)/365      t_exp = (di/d_nom - 1.0) / b / di     qi_exp = qi*np.power(1.0 + b*di*t_exp, -1.0 / b) * np.exp(d_nom*t_exp)      # Separate indices for hyperbolic and exponential decline     mask_hyp = t &lt;= t_exp     mask_exp = ~mask_hyp      # Calculate q_hyp and q_exp for their respective regions     time = np.array(t)     q_hyp = np.zeros_like(time, dtype=float)     q_exp = np.zeros_like(time, dtype=float)      q_hyp[mask_hyp] = arps(qi, di, b, time[mask_hyp])     q_exp[mask_exp] = arps(qi_exp, d_nom, 0.0, time[mask_exp])      # Combine results based on the condition     return q_hyp + q_exp  qi = 1000 di = 0.1 b = 0.5 t = [x for x in range(1, 3650)] lim = 0.5  q_mod_sec = modified_arps(t, qi, di, b, lim, type = 'Secant') q_mod_tan = modified_arps(t, qi, di, b, lim, type = 'Tangent')  plt.plot(t, q_mod_sec, '-', label='Secant Effective') plt.plot(t, q_mod_tan, '-', label='Tangent Effective') plt.legend(loc='upper right') plt.title('DCA') plt.xlabel('Time') plt.ylabel('Production') plt.yscale("log") plt.plot()

В данном листинге приведены примеры реализации модифицированной модели Arps для вычисления кривой производительности на основе эффективного спада по касательной type = ‘Secant’ и эффективного спада по секущей type = ‘Tangent’. В зависимости от типа спада мы вначале вычисляем номинальный наклон. Затем вычисляем период времени, когда необходимо перейти с гиперболической (b > 0) на экспоненциальную (b = 0) модель — t_exp. А также рассчитываем начальную производительность для экcпоненциального периода qi_exp. Расчет производительности можно выполнить отдельно для каждого периода времени, используя индексные маски mask_hyp и mask_exp. Или посчитать q всего периода для b > 0 и b = 0, а затем отфильтровать по времени с помощью выражения np.array(x > t_exp).

t_ext = np.array(x > t_exp) t_exp.astype('int') q = q_hyp*(1 - t_exp) + q_exp*t_exp

В этом примере значение времени (в днях) перехода с гиперболического спада на экспоненциальный равно t_exp = 861.188 — для ‘Secant’ и t_exp = 1033.167 для ‘Tangent’.

Расчет производительности скважины на основе датасета North Dakota Industrial Commission (NDIC)

В этом анализе используется общедоступные данные о добыче из Отдела нефти и газа North Dakota Industrial Commission (NDIC). Файлы с данными можно скачать по ссылке machine-learning-oil-gas-industry. В материалах «Machine Learning in the Oil and Gas Industry» уделяется достаточно внимания подготовке исходных данных (Data preprocessing) для дальнейшего построения модели Arps. В результате получаются два набора данных: train_prod.csv и test_prod.csv. На практике preprocessing зависит от источника данных, например DRI File Spec, и функциональных требований. Это может быть листинг в Jupyter notebook или автоматизированная ETL-система.
В дальнейшем я опускаю шаг обработки данных и использую уже подготовленные наборы train_prod.csv и test_prod.csv.
Загрузим исходные наборы данных и преобразуем ReportDate в формат даты-времени

import pandas as pd  def pre_process(df, column):     df.drop("Unnamed: 0", axis=1, inplace=True)     df.info()     print(df.columns)     # descriptive statistics     df.describe().T     df.head(15)     df.nunique()         df.dtypes     df.shape     # filtering     df.dropna(inplace=True)     # drop rows where oil rate is 0     df = df[(df[column].notnull()) & (df[column] > 0)]     return df  train_prod = pd.read_csv('./train_prod.csv') test_prod = pd.read_csv("./test_prod.csv")  # Basic Processing and data exploration train_prod = pre_process(train_prod, 'Oil') test_prod = pre_process(test_prod, 'Oil')  # convert time to datetime and set as dataframe index train_prod["ReportDate"] = pd.to_datetime(train_prod["ReportDate"]) test_prod["ReportDate"] = pd.to_datetime(test_prod["ReportDate"])

Вычислим производительность нефти и подготовим тренировочный и тестовый набор данных.

df_temp = train_prod[train_prod['API_WELLNO'] == 33025021780000] df_temp["oil_rate"] = df_temp["Oil"] / df_temp["Days"] df_temp['date_delta'] = (df_temp['ReportDate'] - df_temp['ReportDate'].min())/np.timedelta64(1,'D')  df_temp = df_temp[['date_delta', 'oil_rate']] data = df_temp.to_numpy() # T is number of days of production - cumulative # q is production rate T_train, q = data.T  df_temp_test = test_prod[test_prod['API_WELLNO'] == 33025021780000] df_temp_test["oil_rate"] = df_temp_test["Oil"] / df_temp_test["Days"] df_temp_test['date_delta'] = (df_temp_test['ReportDate'] - df_temp_test['ReportDate'].min())  / np.timedelta64(1,'D')  df_temp_test = df_temp_test[['date_delta', 'oil_rate']] data_test = df_temp_test.to_numpy() # T is number of days of production - cumulative # q is production rate T_test, q_test = data_test.T

Чтобы сопоставить исторические данные и выполнить прогноз, используется библиотека с открытым исходным кодом SciPy. Библиотека предоставляет функцию curve_fit(), которая использует нелинейный метод наименьших квадратов для определения необходимых коэффициентов. Метод curve_fit принимает как минимум три параметра $popt, pcov = cerve\_fit(func, T, Q)$ . Первый аргумент в вызове — это анализируемая функция. В нашем случае функцией является Arps модель, которая принимает независимую переменную — время спада, и коэффициенты. Например, в гиперболической функции передаются qi, b и di, тогда как в экспоненциальной функции требуются только qi и di, поскольку b = 0. Норма L2 определяет качество подгонки, где фактическое значение сравнивается с прогнозируемым. Параметр p0 — список начальных значений для подгонки, т. е qi = max(q) (максимальное значение из выборки), di = 0.001, b = 0. Функция curve_fit возвращает две переменные: popt — это наши оптимальные коэффициенты и pcov — расчетная ковариация.

from scipy.optimize import curve_fit  p0 = [max(q), 0.001, 0.0] popt, pcov = curve_fit(arps, T_train, q, p0 = p0) print(popt)

Таким образом, получив коэффициенты, мы может произвести расчет прогнозируемых значений производительности нефти и сравнить их с исходными данными. Где qi = popt[0], di = popt[1], b = popt[2]. В качестве расчетной модели используется modified_arps с эффективным спадом по секущей по умолчанию.

time = pd.date_range(start='6/1/2015', periods= 54, freq='MS') time T_Test2 = T_train[-1] + T_test len(T_train) pred_hyp =  modified_arps(popt[0], popt[1], popt[2], T_train) pred_hyp2 =  modified_arps(popt[0], popt[1], popt[2], T_Test2) print(pred_hyp) print(pred_hyp2) # forecast q_orig = np.append(q, q_test) forecast = np.concatenate([pred_hyp, pred_hyp2])  # hyperbolic forecast - plot plt.plot(time, q_orig, color="black", label='Actual Data') plt.plot(time, forecast, color="green", label="Hyperbolic Trend", linestyle="--") plt.title('Production Forecast') plt.xticks(rotation=45) plt.xlabel('Time') plt.ylabel('Oil Production Rate (bbl/d)') plt.legend(loc='best') plt.show()

Посчитаем среднеквадратическую ошибку и норму L2 с помощью библиотеки sklearn.

from sklearn.metrics import r2_score from sklearn.metrics import mean_squared_error from math import sqrt  rmse = sqrt(mean_squared_error(q_orig, forecast)) print("RMSE - Hyperbolic Method:", rmse) #RMSE - Hyperbolic Method: 81.40806677104752  r2_score(q_orig, forecast) # 0.9496932527194057

Подобным образом может выполняться прогнозирование производительности скважины, используя альтернативные DCA модели. Где в качестве модельной функции для подгонки параметров могут служить методы: SEPT, Duong, PowerLow.

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

ETL (extract load transformation) — система. Принимает исходные данные из внешних систем, выполняет очистку и подготовку данных, загружает обработанные данные в базу
Сервис с алгоритмами прогнозирования. Может автоматически анализировать новые данные и выполнять обработку в многопоточном режиме.
База данных с обработанными входными данными и вычисленным прогнозом
Frontend-сервис, где пользователь получает визуальную информацию по историческим и прогнозируемым данным

Применение алгоритмов машинного обучения для прогнозирования производительности добычи нефти из скважин позволяет решать множество актуальных задач. Использование современных программных инструментов для разработки моделей машинного обучения предоставляет широкие возможности благодаря богатой экосистеме библиотек, таких как Scikit-learn, SciPy, Pandas, Numpy и других. Эти технологии позволяют эффективно анализировать большие объемы данных, выявлять скрытые закономерности и оптимизировать производственные процессы.

Используемые источники:

ссылка на оригинал статьи https://habr.com/ru/articles/865442/