CC BY
43
Известия Тульского государственного университета Естественные науки. 2009. Вып. 3. С. 255-266
= Науки о земле =
УДК 622.276.5.001.42+004.89
Решение обратных коэффициентных задач с применением нечетких деревьев регрессии на примере обработки кривой восстановления давления
Р.Н. Исмагилов, В.Е. Лялин, К.А. Сидельников
Аннотация. Авторами предложено для обработки результатов гидродинамических исследований скважин (ГДИС) использовать систему нечеткого логического вывода по Сугено, база правил которой представляет собой нечеткое дерево регрессии. Сформулированы основные этапы построения таких нечетких систем. Приведены результаты обработки кривой восстановления давления.
Ключевые слова-, гидродинамические исследования скважин, идентификация системы, деревья классификации и регрессии, нечеткая система.
В течение ГДИС осуществляется мониторинг реагирования пласта на изменение условий добычи (нагнетания). Поскольку эта реакция в какой-то степени характеризует свойства пласта, то последние во многих случаях можно, так или иначе, оценить. Таким образом, интерпретация результатов ГДИС представляет собой обратную задачу, при которой параметры модели определяются путем анализа отклика модели на входное воздействие.
Математические модели пласта-коллектора создаются путем решения диффузионного уравнения при различных граничных условиях. Модели определяются тремя различными компонентами, которые характеризуют пласт, скважину и ее окрестности (внутренние граничные условия) и внешние границы пласта (внешние граничные условия). При определенных допущениях фильтрация слабо-сжимаемого флюида в пласте подчиняется диффузионному уравнению
д2Р 1 дР фцсгдР т
дг2 + г дг ~ к '
где Р — давление; г — расстояние по радиусу от ствола скважины; I —
время; ф — пористость; /х — вязкость; с% — общая сжимаемость системы; к —
абсолютная проницаемость.
Как правило, аналитическое решение уравнения (1) проще (а иногда и единственно возможно) получить в пространстве изображений, чем во временной области. Например, для модели бесконечного пласта, характеризуемой тремя параметрами к, Б и С, изменение безразмерного забойного давления можно определить по следующей формуле [1]:
где Р£> — безразмерное давление в пространстве изображений; Со — безразмерный коэффициент накопления; Ко, К\ — модифицированные функции Бесселя второго рода 0-го и 1-го порядков соответственно.
Для оценки параметров во временной области требуется обратное преобразование Лапласа. В частности, обратное преобразование уравнения (2) должно быть выполнено численными способами. В [2] описан ряд алгоритмов для решения этой задачи: Гавера-Стехфеста, Эйлера и Талбота. Последние два метода, хотя и быстрее первого, однако требуют операций с комплексными числами. Подобные алгоритмы, в общем, базируются на вычислении следующей линейной комбинации значений изображения [2]:
где / — неотрицательная вещественная функция (оригинал); / — изображение /; и ск& — соответственно вес и узел, которые представляют собой комплексные числа, зависящие от п, но не зависящие от / или времени Расчеты по формуле (3) требуют достаточно высокой машинной точности.
Преобразование Лапласа является интегральным преобразованием, что приводит к сглаживанию погрешностей экспериментальных функций. Кроме того, часто, исходя из точного решения в изображениях, удается получить асимптотики при £ —>• оо и £ —>• 0 решений во временной области и эффективно использовать их при решении обратных задач [3].
Так или иначе, применение преобразования Лапласа ограничено, очевидно, только теми случаями, когда решение (1) может быть получено аналитически. Однако это возможно далеко не всегда. Поэтому если точное решение получить не удается, можно произвести замену дифференциального уравнения (1) его конечно-разностным аналогом [3]. Определение параметров линейных и нелинейных дифференциальных уравнений относится к построению так называемой модели «серого ящика». Модель называется так, поскольку ее математическая структура задается в явном виде. При этом взаимосвязи между всеми параметрами и переменными модели известны. Модель серого ящика удобна для использования различных методов идентификации систем [4].
Рассмотрим пример простой модели распространения тепла в изолированном металлическом стержне длиной Ь. Один конец стержня нагревается
Рв
Ко (у/з) + Бу/зКг (у/з)
(2)
я {у/зКг (у/з) + Сов {Ко (у/з) + ву/зКг (у/в))) ’
О < * < оо,
(3)
&=о
с тепловой мощностью и (¿), на другом конце осуществляется измерение температуры у (¿). При идеальных условиях данная система описывается уравнением теплопроводности (диффузии)
дх (і, С) <92ж(і,£) — к-
ді
де
(4)
где к — коэффициент теплопроводности; ж(£,£) — температура стержня в точке £ в момент времени
Для получения непрерывной во времени модели в пространстве состояний вторая производная в правой части (4) заменяется конечно-разностной аппроксимацией
в2х (*, С) ж (г, £ + АЬ) - 2ж (*, £) + ж (г, £ - АЬ)
(АЬУ
(5)
где £ = кАЬ; А: = 1,2,... ,п; п = Ь/АЬ. Подобное преобразование приводит к системе с конечным набором п состояний ж (і, кАЬ), каждое из которых представляет собой некую усредненную температуру стержня в диапазоне кАЬ ^ £ < (к + 1) АЬ. Другими словами порядок системы зависит от размера АЬ пространственной сетки.
В результате, уравнение диффузии (4) можно представить как непрерывную во времени линейную систему в пространстве состояний:
х(і) = Ах (і) + Ви (і) + Ке (і), х(0) = хо, у (і) = Сх (і) + Би (і) + е(і),
(6)
где х (£) — вектор состояний; и (¿) — вектор входных воздействий; е (¿) — вектор шумов; у (£) — вектор реакций системы; хо — вектор начальных состояний системы. Для нашей задачи
уи
ж (і, 0) ж (0,0)
х ж (і, кАЬ) х ж (0, кАЬ)
і * * ч-Г н 1 х (0, Ь)
А
АЬ2
'-1 1 0 . . 0 0 ' '1 /АЬ' '0'
1 -2 1 . . 0 0 0 0
0 1 -2 . . 0 0 0 0
В К
0 0 0 . . -2 1 0 0
. 0 0 0 . . 1 -1. 0 0.
С = [1 0 0 ... О 0] ; Б = [0]
Далее, имея достаточный объем наблюдений, можно использовать любой подходящий метод идентификации систем (например, метод минимизации ошибки прогноза (РЕМ), который во многом основан на авторегрессионной модели скользящего среднего с учетом внешних воздействий (ARMАХ)) для оценки параметра к.
Несмотря на большую гибкость описанного способа замены дифференциального уравнения конечно-разностным, у него есть ряд серьезных недостатков. Данный подход следует применять с большой осторожностью, поскольку конечно-разностная аппроксимация эквивалентна непосредственному дифференцированию экспериментальных функций, чреватому большими погрешностями [3]. Кроме того, существует, например, модель пласта с непроводящим сбросом, которая, по сути, является двумерной. Если использовать преобразование Лапласа, то применительно к данной модели существует аналитическое выражение для изменения давления в скважине [1]. Однако использование численного решения этой задачи для получения системы (6) может потребовать чрезмерно большого количества переменных состояния. Также можно отметить, что существуют определенные трудности в создании модели вида (6) для обработки данных ГДИС, полученные методом восстановления давления. Проблема связана с зависимостью вектора начальных условий от параметров пласта.
Помимо упомянутых трудностей, возникающих при работе с конечными разностями, существует, пожалуй, главное препятствие для более широко применения моделей типа (6). Для каждого нового набора наблюдений требуется каждый раз численно решать соответствующее дифференциальное уравнение, что связано с большими вычислительными затратами. Удобней было бы заранее иметь некоторое множество численных решений, по которым можно было бы аппроксимировать сложную функциональную зависимость между независимыми переменными (в данном случае расход скважины), параметрами (в данном случае пласта) и зависимыми переменными (в данном случае давление скважины). Для подобных целей широко используются так называемые методы мягких вычислений (soft computing): нейронные сети и системы нечеткого логического вывода. Как показывает ряд исследований [5, 6], последние обычно предпочтительнее. Причинами являются более простая настройка нечетких систем, а также более понятная внутренняя структура и логика работы таких систем, основанная на правилах.
Современные системы нечеткого логического вывода основываются на применении того или иного метода структурной оптимизации для определения структуры ее базы правил. К таким методам относятся различные способы кластеризации или классификации исходных данных, выступающих в роли обучающей выборки для нечеткой системы. Так, например, достаточно успешно применяется нечеткая кластеризация по методу С-средних (fuzzy С-means или FCM). Тем не менее, для решения задач нелинейной регрессии особую популярность приобрел метод деревьев классификации и регрессии (classification and regression tree или CART), являющимися частными
случаями деревьев решений. С небольшими изменениями CART можно использовать для идентификации структуры базы правил нечеткой системы. Рассмотрим основные шаги, необходимые для генерации системы нечеткого логического вывода на базе CART применительно к анализу данных ГДИС.
1. Формирование обучающей и тестовой выборок
Данные ГДИС методом восстановления давления формировались на основе аналитической модели (2):
Pdbu (*£>) = Pd (tpD) + Pd (*£>) — Pd (tpD + to), (7)
где Pdbu — безразмерное изменение давления при ГДИС методом восстановления давления; I’d — безразмерное изменение давления при ГДИС методом понижения уровня (определяется по (2); to — безразмерное изменение времени; tpD — безразмерное время, в течение которого осуществлялась добыча до закрытия скважины.
Информация о пласте и насыщающем его флюиде взята из [7]. Общее число точек Nt за период времени от 0,01 до 100 часов составило 51. Всего случайным образом генерировалось N,. = 50 кривых изменения давления. Половина из них использовалась в качестве обучающей выборки, другая половина — в качестве тестовой. При этом областями изменений параметров k, S и С являлись соответственно [0, 01; 0,1] мкм2, [0; 10] и [0,01; 0,1] м3/ МПа. Итого, суммарный объем точек наблюдений Ns = NtNc = 2550.
Ввиду большого разброса значений времени и изменения давления, вместо реальных величин t и А Рви ПРИ обучении использовались соответственно log (t) и log (АРви)- Это положительно повлияло на весь процесс обучения, так известно, что правильное масштабирование зависимых и независимых переменных может существенно улучшить обусловленность задачи [8].
2. Построение оптимального дерева регрессии
В соответствии с алгоритмом CART строилось дерево регрессии для прогноза отклика А Рви как функции предикторов к, S, С и t. К задаваемым пользователем параметрам алгоритма относятся:
— целое число minparent, которое устанавливает для каждого узла дерева минимальное количество наблюдений, участвующих в разбиении (равнялось 10);
— целое число minleaf, которое устанавливает для каждого листа дерева минимальное количество содержащихся в нем наблюдений (равнялось 1);
— флаг mergeleaves, который разрешает или запрещает объединять листья, имеющие одного и того же родителя, если сумма соответствующих им рисков равна или превышает значение риска для родительского узла (разрешал);
— вещественное число де£о/ег, которое определяет для каждого узла минимальный относительный порог квадрата ошибок, влияющий на решение продолжать или нет процесс разбиения в данном узле (равнялось 10~6).
20 40 60 80 100 120 140 160 180 200
Размер дерева (число листьев)
Рис. 1. Стоимость дерева в зависимости от его размера
Чтобы избежать переобучения дерева, проводилось его тестирование с выбором оптимального уровня отсечения. Для этого рассчитывалась стоимость дерева как взвешенная сумма стоимостей узлов. Весом выступает вероятность попадания в данный узел. Стоимость узла вычисляется как среднеквадратическая ошибка по всем наблюдениям, находящимся в этом узле.
В результате был построен график на рис. 1, где пунктирная линия соответствует значению стоимости дерева, большей минимальной (помечена кружком) на одну величину среднеквадратической ошибки. Там же квадратом. помечено дерево (уровень отсечения) с наименьшим количеством листьев и стоимостью, расположенной ниже пунктирной линии. Это дерево с 35 листьями и было выбрано для дальнейших расчетов.
3. Генерация базы правил системы нечеткого логического вывода
В работе нечеткая система осуществляла логический вывод по Сугено (в частности, используется нечеткая модель Сугено первого порядка), так как это упрощает дальнейшую параметрическую настройку системы.
Поскольку дерево регрессии — это частный случай дерева решений, то его можно представить в виде набора четких правил. Обозначим множество независимых переменных и параметров на входе дерева как { | г = 1, /}
и зависимую переменную на выходе дерева как у. Любое правило состоит из посылки и следствия. Посылка в данном случае представляет собой результат сравнения Жг с некоторым пороговым значением тп в п-м узле дерева (п = 1, М). Пусть г?п Е 1,/ — это индекс входной переменной, по которой
происходит сравнение в п-м узле дерева. У каждого узла дерева, кроме терминальных узлов (листьев), есть два потомка: левый и правый. Если результат ж„п < гип истинен, то дальнейшие сравнения проводятся в левом дочернем узле. Если истинен результат ж„п ^ пзп, то дальнейшие сравнения проводятся в правом дочернем узле. В дальнейшем, не конкретизируя тип отношения, будем записывать просто ж„п о иіп. Тогда любое правило можно представить в следующем общем виде:
где т 6 1, N — индекс узла, только через который можно попасть в 1-й лист.
Для того чтобы четкое правило (8) сделать нечетким, необходимо фаз-зифицировать посылки и следствие правила. В посылке ж„п о щп заменяется операцией принадлежности нечеткому множеству Лп, т. е. ж„п 6 Лп■ Это равносильно вычислению функции принадлежности [лдп (ж„п), в роли которой была выбрана сигмоидная функция:
где обычно сп = гип; ап либо настраивается пользователем, либо вычисляется автоматически. Причем знак ап определяет тип нечеткого отношения.
В процессе построения дерева использовалась обучающая выборка, которую можно представить как упорядоченное множество наблюдений =
числе и терминальному, можно установить подмножество С 5 как результат разбиения С 5 для п-го узла, являющегося родителем с/-го узла. Разбиение осуществляется в соответствии с ж„п о щп. При этом справедливы равенства
В результате параметр ап в (9) можно связать с ап следующим образом:
ЕСЛИ
(8)
т
)
)
і
(9)
Я4п (ж; ап, Сп) = вщ (ж; ап, сп) =
1 -)- е~ап(х-сп) ’
{ (®м, • • • , Ж/,ё, Уа)I (1 = 1, £>}. При обходе дерева каждому д-му узлу, в том
Вычислим следующую среднеквадратическую ошибку:
(10)
(П)
где ег£(ж) — функция ошибок; К > 0.
Выбор ап по (11) обеспечивает равенство площадей под sig (ж; ап, сп) и функцией Лапласа Ф (ж; сп, сгп) в интервале [— Кап; Кап]. Эксперименты показали, что выбор К = 1 вполне удовлетворителен.
Следствие в нечетком правиле по Сугено представляется в следующем виде:
/
fi=Po,i + ^2pi,iXi. (12)
i—1
Обозначим X; = [ж^]Т — матрица наблюдений входных переменных, попавших в /-й лист, и yi = [у^] - вектор наблюдений выходной переменной в этом же листе. Тогда для нахождения вектора неизвестных параметров РI = [Pm,l], т = О, I требуется решить следующие нормальные уравнения:
XfXm = Xfyh (13)
где Xr = [X| 1] — дополненная единицами справа матрица входных наблюдений.
Решение (13) может осуществляться любым подходящим методом фак-XX воспользоваться моделью Сугено нулевого порядка:
Ро,1 = У1, Pi,l = 0, i = 1,1. (14)
Следует сказать несколько слов о структуре нечетких правил, получаемых на базе CART. Число посылок каждого правила может превышать число входных переменных. Кроме того, любая входная переменная может входить в правило несколько раз. То есть структура нечетких правил на базе CART или любого другого дерева решений достаточно сложная и неоднородная. Однако основная проблема даже не в этом. В наиболее популярных программных пакетах, реализующих системы нечеткого логического вывода, таких как FuzzyTECH или Fuzzy Logic Toolbox в MATLAB принципиально невозможно реализовать правила вида (8). Это связано с негибкой формой задания базы правил в этих пакетах. Поэтому авторам, использовавшим MATLAB для получения приведенных здесь результатов, пришлось модифицировать стандартные модули (М-файлы, листинги на языке С) для работы с подобными правилами.
С другой стороны, если бы мы генерировали базу правил с использованием FCM, то в этом случае правила имели бы стандартный вид. Однако при одинаковом количестве правил, генерируемых CART и FCM, у первых общее количество нелинейных параметров (то есть параметров функций принадлежности) обычно в разы меньше, чем у вторых. Объясняется это тем, что одна и та же функция принадлежности некоторого входа может входить в разные правила, отличаясь только знаком при а. Данное обстоятельство может существенно упростить выполнение следующего этапа.
4. Обучение системы нечеткого логического вывода
После того как будет определена структура нечеткой системы требуется ее параметрическая настройка (оптимизация). Можно выделить два подхода к решению данной задачи как частного случая нечеткой идентификации (fuzzy identification) параметров.
В первом случае оптимизируются только линейные параметры нечеткой базы правил — коэффициенты в заключениях правил [9]. Для этого можно использовать псевдоинверсию Мура-Пенроуза, фильтр Калмана (рекуррентный метод наименьших квадратов) или градиентный метод настройки. Чаще используется фильтр Калмана, так как он дает более устойчивые оценки параметров по сравнению с псевдоинверсией и не требует задания начальных приближений параметров в отличие от градиентного метода. Тем не менее встречаются исследования [10], в которых авторы успешно применяют псевдоинверсию вкупе с тем или иным методом регуляризации.
Число эпох
Рис. 2. Среднеквадратическая ошибка по обучающей и тестовой выборкам
Во втором случае одновременно настраиваются и заключения правил, и функции принадлежности термов входных переменных. Особенно хорошо себя зарекомендовал гибридный метод оптимизации ANFIS [6] (адаптивная нейро-нечеткая системы логического вывода), в котором линейные параметры настраиваются по фильтру Калмана, а нелинейные параметры — по градиентному методу обратного распространения ошибки. Однако, в принципе, все параметры можно настраивать только по градиентному методу, то есть подобно алгоритмам обучения нейронных сетей.
В работе настройка осуществлялась на базе ANFIS с использованием той же обучающей выборки, что и при построении CART. Во избежание переобучения нечеткой системы на данном этапе проводился контроль по тестовой выборке. На рис. 2 приведены графики среднеквадратических ошибок, полученных по обучающей и тестовой выборкам. Квадратом отмечена нечеткая
система, имеющая минимальную ошибку на тестовом множестве. На рис. 3 построен кросс-плот, где по оси абсцисс расположены нормализованные значения, выдаваемые нечеткой системой, а по оси ординат — нормализованные значения из тестовой выборки. Коэффициент корреляции р = 0,9968.
Рис. 3. Кросс-плот нормализованных фактических и расчетных значений
Оценки параметров и дисперсии
Модель П роницаемость &, мкм2 Скин-фактор S Коэффициент накопления С, м3/МПа Дисперсия <т2, МПа2
Аналитическая модель 0,0499 ±0,0146 % 4,9819 ±0,3569 0,0500 ±0,0023 % 0,8848- Ю”'4
Нечеткое CART 0,0492 ±0,0150 % 4,9307 ±0,3096 0,0506 ±0,0064 % 23,71 ■ 10“1
5. Применение нечеткой системы для обработки данных ГДИС
Для имитации случайных ошибок, распределенных по нормальному закону, в данных ГДИС использовался генератор случайных чисел с нулевым математическим ожиданием и дисперсией 10-4 МПа2. Истинные значения параметров были к = 0,05 мкм2, S = 5, С = 0,05 м3/МПа.
Результаты оценок параметров с помощью аналитической модели и нечеткого CART приведены в таблице. Полученные обоими методами кривые изменения давления и его производной представлены на рис. 4.
На основании проведенных исследований можно сделать следующие выводы. Нечеткие системы, структура базы правил которых формируются на базе CART, действительно способны аппроксимировать сложные нелинейные зависимости между независимыми переменными, параметрами и зависимыми
Рис. 4. Динамика изменения давления и его производной
переменными. Подобные системы достаточно успешно можно использовать для проведения нелинейного регрессионного анализа. Однако от них не следует ожидать высокой точности получаемых оценок, как в случае аналитических или численных моделей. Основное применение систем нечеткого логического вывода видится в проведение различий между моделями-кандидатами. Например, нечеткие модели пласта, обученные на численных решениях задач притока нефти или газа к скважине, могут стать основой дискриминантного анализа по методу последовательного прогнозирования вероятностей [7]. Выбрав подходящую модель пласта, оценки параметров затем можно уточнить с помощью любых методов, описанных в начале статьи.
Список литературы
1. Anraku Т. Discrimination between reservoir models in well test analysis. Ph.D. Thesis. Stanford University, 1993.
2. Abate JWhitt W. A unified framework for numerically inverting Laplace transforms // INFORMS J. on Computing. 2006.
3. Мирзаджапзаде A.X., Хасанов M.M., Бахтизин P.H. Моделирование процессов нефтегазодобычи. Нелинейность, неравновесность, неопределенность. М.-Ижевск: Институт компьютерных исследований, 2005. 368 с.
4. Гроп Д. Методы идентификации систем. М.: Мир, 1979. 302 с.
5. Тененев В,А., Ворончак В,И. Решение задач классификации и аппроксимации с применением нечетких деревьев решений // Интеллектуальные системы в производстве. Ижевск: ИжГТУ, 2005. №2. С. 46-69.
6. Jang J.-S.R., Sun С.-Т., Mizutani Е. Neuro-fuzzy and soft computing. Prentice-Hall, 1997.
7. Сидельпиков К.А., Денисов C.B. Метод последовательного прогнозирования вероятностей для дискриминатного анализа моделей пласта по данным ГДИС // Искусственный интеллект. Интеллектуальные системы: матер. IX Межд. научно-техн. конф. / Донецк: ИПИИ «Наука i ocBÎTa», 2008. T. 2. С.74-83.
8. Дэннис Док., Шнабель Р. Численные методы безусловной оптимизации и решения нелинейных уравнений. М.: Мир, 1988. 440 с.
9. Штовба С.Д. Проектирование нечетких систем средствами MATLAB. М.: Горячая линия — Телеком, 2007. 288 с.
10. Johansen Т.A. Robust identification of Takagi-Sugeno-Kang fuzzy models using regularization // Fifth IEEE International Conference on Fuzzy Systems. New Orleans, 1996.
Поступило 05.09.2009
Исмагилов Рустам Наилевич (gk34@yandex.ru), начальник технического отдела, Уренгойское газопромысловое управление, ООО «Газпром добыча Уренгой».
Лялин Вадим Евгеньевич (velyalin@mail.ru), д.т.н., профессор, заведующий кафедрой, кафедра интеллектуальных информационных технологий, Ижевский государственный технический университет.
Сидельников Константин Анатольевич (sidelkin@yandex.ru), к.т.н., доцент, кафедра интеллектуальных информационных технологий, Ижевский государственный технический университет.
Parameter estimation using fuzzy regression tree in terms of pressure build-up curve processing
R.N. Ismagilov, V.E. Lyalin, K.A. Sidelnikov
Abstract. Authors propose Sugeno fuzzy inference system with rule base, represented by fuzzy regression tree, to process well test data. The main steps of such fuzzy systems construction are described. The results of pressure build-up curve processing are presented.
Keywords: well test, system identification, CART, fuzzy system.
Ismagilov Rustam (gk34@yandex.ru), Gazprom dobycha Urengoi Ltd.
Lyalin Vadim (velyalin@mail.ru), doctor of technical sciences, professor, chair of department, department of Intelligent IT, Izhevsk State Technical University.
Sidelnikov Konstantin (sidelkin@yandex.ru), candidate of technical sciences, associate professor, department of Intelligent IT, Izhevsk State Technical University.
