ЗАДАЧА СТРУКТУРНО-ПАРАМЕТРИЧЕСКОЙ ИДЕНТИФИКАЦИИ СИСТЕМ ПРИ МОДЕЛИРОВАНИИ ДВУХФАЗНОЙ ФИЛЬТРАЦИИ В ПОРИСТЫХ СРЕДАХ Текст научной статьи по специальности «Математика»

УДК 622.276.031:531.72.5:519.87

ЗАДАЧА СТРУКТУРНО-ПАРАМЕТРИЧЕСКОЙ ИДЕНТИФИКАЦИИ СИСТЕМ ПРИ МОДЕЛИРОВАНИИ ДВУХФАЗНОЙ ФИЛЬТРАЦИИ В ПОРИСТЫХ СРЕДАХ

Э. Н. Мусакаев 123, С. П. Родионов 123, В. П. Косяков 123

1 Тюменский филиал института теоретической и прикладной механики им. С. А. Христиановича Сибирского отделения Российской академии наук

2 ООО «Юни-КОНКОРД»

3 Тюменский государственный университет, musakaev91@gmail.com, rodionovsp@bk.ru, hammer-rav@mail.ru

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

Ключевые слова: модель, фильтрация, пористая среда, структурно-параметрическая идентификация.

PROBLEM OF STRUCTURE AND PARAMETER IDENTIFICATION AT SIMULATION OF TWO-PHASE FILTRATION IN POROUS MEDIA

E. N. Musakaev 123, S. P. Rodionov 123, V. P. Kosyakov 123

1 Khristianovich Institute of Theoretical and Applied Mechanics, Tyumen site, Siberian Branch,

Russian Academy of Sciences 2 OOO Uni-CONCORD

3 Tyumen State University, musakaev91@gmail.com, rodionovsp@bk.ru, hammer-rav@mail.ru

The article proposes a technique for the structure and parameter identification of systems described by models of single-phase and two-phase filtration in a porous medium with the presence of strong heterogeneities. There are some formulated principles that need to be followed during constructing the model structure. The paper also includes an algorithm for detecting fractures and sealing faults using only production data. A good agreement of the results obtained on the reference model and on the model constructed with the use of the proposed structure and parameter identification technique is shown.

Keywords: model, filtration, porous medium, structural-parameter identification.

Введение

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

мы дискретизированных конечно-разностных уравнений фильтрации, приводящая к снижению качества прогнозных величин, таких как дебит жидкости и дебит нефти. Если численная и/или математическая модель содержит слишком много неизвестных параметров (т. е. является слишком сложной), то, после выполнения параметрической адаптации, такая модель не обладает устойчивостью прогноза [1, 2]. Поэтому актуальной становится задача структурной идентификации модели фильтрации, т. е. задача выбора фильтрационной модели оптимальной сложности, которая, с одной стороны, описывала бы процесс достаточно точно, а с другой - обладала бы хорошими прогнозными свойствами.

В большинстве работ по синтезу математических моделей их структура, как правило, постулируется (выбирается) априори с точностью до некоторого множества неизвестных параметров. В дальнейшем это множество является основным объектом исследования. Доминирующим подходом при выборе структуры модели в настоящее время является статистический подход, когда решение принимается на заданном классе моделей-претендентов [3].

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

Методика структурно-параметрической идентификации

Фильтрационная модель нефтяного пласта обычно содержит большое количество расчетных ячеек. Поэтому заполнение необходимых массивов данных о свойствах породы и жидкостей требует, во-первых, решать проблему интерполяции и экстраполяции данных измерений на скважинах на межскважинное пространство, а во-вторых, решать проблему усреднения или масштабирования данных, полученных на масштабах керна и геофизических исследований, на масштаб расчетных блоков [4].

В условиях высокой неопределенности исходных данных и большой погрешности измерений (особенно это касается межскважинного пространства) нарушается принцип «простоты», который в реалиях математического моделирования сложных объектов может звучать так: «Модель должна быть настолько простой насколько возможно и, настолько сложной насколько необходимо» [5].

Interwell model — фильтрационная модель потоков жидкости между скважинами. Представим гидродинамическую модель нефтяного месторождения, разрабатываемого нагнетательными и добывающими скважинами, как модель с парными связями между скважинами (рис. 1). Как показано на рис. 1а красными линиями, скважины (или любые объекты с заданными границами, например, участки) имеют «объем дренирования» Vt. Фильтрация флюидов между скважинами с номерами i и j характеризуется двумя параметрами: проводимостью Tij и объемом дренирования Vp,ij, внутри которого проводимость определяет потоки флюидов. Определим проводимости между скважинами следующим образом:

_ BijKjjHjj li,j = Liji4j '

где К - абсолютная проницаемость на границе между скважинами; В - ширина стенки, через которую фильтруется жидкость; Н - среднеарифметическая эффективная толщина коллектора; L - расстояние между скважинами; ц - вязкость жидкости.

Вязкость определяется через начальную водонасыщенность:

М- + (1 Smw)№o,

где - начальная водонасыщенность; Мw - вязкость воды; М0 - вязкость нефти.

а

Ь

Рис. 1. Схема построения связей между скважинами и дренируемых объемов:

а - построение дренируемого объема ^ для пары скважин «/-/'»; Ь - для множества связанных скважин

Система уравнений материального баланса. Для идентификации потоков жидкости между скважинами с использованием значений дебитов жидкости и забойных давлений на скважинах применяется система уравнений материального баланса для однородной жидкости [6], которая имеет следующий вид:

№ ^ = кдзЛ - джЛ + ти}(р1 - р}) - № - Рпвд.

(2)

Здесь и далее индексом / внизу (/ = 1, ...,М) отмечены величины, относящиеся к -ой скважине: Р^ - среднее пластовое давление; Qж,i и Qз,I - соответственно объемы добытой жидкости и закачанной воды; - эффективный коэффициент сжимаемости смеси жидкостей и горной породы; к^ - коэффициент полезной закачки воды, равный доле закачанной воды, которая поступает непосредственно в расчетный пласт; Хпв,1 - коэффициент, определяющий интенсивность притока подошвенной воды; щ - количество скважин, связанных с I-ой скважиной; через Рпв,г обозначено среднее давление в области подошвенной воды.

Далее проводится дискретизация системы уравнений (2) и записывается в матричном

виде:

№ ^ = к£зЛ - + ти(Р1 - р)) - хивЛ(Р1 - равЛ), ~ лрп+1 = ь,

где & = tn -^, (п = 1,..., N) - шаг по времени;

(3)

п - верхний индекс, обозначающий переменную в момент времени tn; А - симметричная матрица размерности N X N с диагональными элементами (А)ц =

в^ + Ъ^Л и недиагональными (А)и =-Ти; Ъ = ^Рп + Qn3^1 - Qnжл1 - вектор-

столбец размерности N.

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

Параметрическая идентификация. Задача идентификации параметров модели по промысловым данным относится к классу обратных задач. Если при решении прямой задачи

(3) определяются пластовые давления при заданных параметрах пласта (ß,k,XnB,T), то в обратных задачах находится причинная составляющая наблюдаемых событий с использованием некоторых косвенных данных или априорных знаний о структуре пласта.

Задача идентификации параметров модели к промысловым данным ставится и решается в оптимизационной постановке. Решением задачи будет являться такое значение вектора управляющих параметров u = uovt, при котором целевая функцияJ(u) имеет глобальный минимум j(uopt) ^ min.

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

и = щ,, (4)

Ttj = ацТи, (5)

где aij Е [0, атах]. Если а^ = 0, то связь между скважинами i и j (перетоков жидкости) отсутствует.

Целевую функцию представим в виде меры отклонения фактических и расчетных значений давления в течение истории разработки:

, Гп = 1|^ф>акт ^модель (u)|

Jp(и) =--, (6)

^п=1^факт

где I - число временных шагов;

Й - вектор управляющих параметров;

Рфакт(п = 1, ...,1) - значения пластового давления, рассчитанное с использованием данных промысловой истории в момент времени t = tn;

РМоДелъ(и) - аналогичные значения пластового давления, определяемые из решения уравнений (2) при значении вектора управляющих параметров Й.

Путем идентификации множителей на проводимости (4), получаются адаптированные значения интенсивности потоков жидкости между скважинами (5), где значения (1) рассчитываются на основе входных данных.

Количество идентифицируемых параметров зависит от матрицы связей, определяемой вручную или генерируемой автоматически.

Структурная идентификация и регуляризирующие функционалы. При традиционном подходе к моделированию нефтяного месторождения с помощью коммерческих гидродинамических симуляторов число параметров, которые можно определить из имеющейся геолого-промысловой информации, много меньше количества параметров пласта, необходимых для расчета. Поэтому решение задачи идентификации не уникально или, другими словами, существует не единственный набор параметров, при котором целевой функционал имеет минимум. Усложнение модели путем добавления в нее новых параметров, учитывающих новые особенности реального объекта, может привести приводит к желаемой минимизации разницы между наблюдаемыми и расчетными данными. С другой стороны, это приводит и к проблеме не уникальности при решении задачи параметрической идентификации. Прогноз с помощью такой численной и/или математической модели после параметрической адаптации будет неустойчив или вовсе неверным [5].

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

Решением задачи структурной идентификации будет являться такая модель оптимальной сложности Mopt Е МсШ55, которая подходит под следующие условия:

1. Модель удовлетворительно описывает промысловые данные (качественная адаптация на историю разработки).

2. Модель дает устойчивый прогноз.

Первое условие выполняется в большинстве случаев благодаря параметрической идентификации, для выполнения второго условия важно учесть следующие пункты:

1. Необходимо выделять обучающую и независимую контрольную выборки данных для отбора наилучших моделей.

2. Необходимо выбирать модель оптимальной сложности.

Рассмотрим каждый пункт более подробно и более конкретно к рассматриваемой задаче.

1. Период разработки месторождения Та11 разделяется на период «обучения» То и период «экзамена» Те,

Та11 =То +Те.

На интервалах, относящихся к периоду «обучения», производится параметрическая идентификация и, при помощи формулы (6), будет получено значение функционала J0(u0pt). Теперь, зная 1 = uopt и решая прямую задачу (3) на периоде «экзамена», получим значение функционала Jе. Общий функционал будет тогда равен:

Jal l =Jo+ YJe, (7)

где у - весовой коэффициент на период «экзамена».

2. Для учета сложности выбранной модели используем байесовский информационный критерий (BIC):

SSF

BIC = L^n^ + k^nL, (8)

где SSE - сумма регрессионных остатков; L - количество наблюдаемых величин; к - параметр сложности модели.

Формирование класса моделей-претендентов Mclass(lapriori), из множества которых будет производиться выбор модели оптимальной сложности M0pt, обычно происходит за счет априорной информации Iapri0ri, общих соображений специалиста; и здесь ключевую роль начинает играть опыт, интуиция человека, что позволяет значительно сузить круг физико-математических моделей до приемлемого, и рассматривать ограниченный набор структур. По мере решения задачи идентификации, подгонки различных параметров модели, накапливается опыт по данному вопросу, производится сравнение различных структур и растет понимание «истинной» системы [7], как она функционирует, вручную выстраиваются «связи» моделей, например, при описании сложного реального объекта сначала простыми моделями, потом более сложными.

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

Приведем общую блок-схему структурно-параметрической идентификации (рис. 2):

Задание начального множества моделей-претендентов Mclass (Iapriori)

Блок структурной идентификации

Формирование вектора управлений и для текущей модели M Е Mclass

Блок параметрической идентификации

Минимизация функционала на период обученияJ 0 = J0(uopt) ^ min

Решение прямой задачи на период экзамена Je =J e(uoVt)

Учет сложности выбранной модели BIC = BIC(uopt,M)

Синтез функционала параметрической идентификации ] p = Jp(BI C,J0,Je)

В соответствие текущей модели M ставится функционал ] p

Минимизация функционала Jstr = J str(Mopt'Jp)

Блок статистического анализа

Выборка лучших Nm моделей с точки зрения функционала ] s tr

Статистическая обработка этих моделей и получение новой информации Istat об объекте моделирования

Сужение множества моделей-претендентов Mclass Mclass(Japriori,^stat)

Рис. 2. Блок-схема метода структурно-параметрической идентификации

Пример применения методики структурно-параметрической идентификации

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

Фильтрационная модель, построенная для расчетов на гидродинамическом симуляторе, будет считаться исходной или «истинной»; качество структурно-параметрической идентификации на настраиваемой модели «Interwell model» будет оценено с помощью «истинной» модели. Для «истинной» модели будет известен весь необходимый массив данных о

свойствах породы и жидкостей, настраиваемая модель в качестве входных данных будет использовать лишь данные по скважинам: промысловые данные и данные РИГИС. Задачей будет идентификация связей между скважинами на настраиваемой модели с применением методики параметрической и структурной идентификации и, таким образом, определение расположения трещины и экрана в исходной модели.

Постановка задачи. Имеется нефтенасыщенный участок (рис. 3-4), который разрабатывается 9 скважинами (5 нагнетательных и 4 добывающие). На участке имеется трещина, расположенная от нагнетательной скважины № 1 до области добывающей скважины № 6. Также имеется протяженная непроницаемая зона, «экран». Параметры задачи указаны в табл. 1.

Рис. 3. Геометрическое расположение неоднородности коллектора на участке

Рис. 4. Распределение пластового давления на участке

Таблица 1

Параметры

Параметр Значение

Размер участка 600х600 м

Пористость 0.3

Толщина пласта 10 (м)

Вязкость воды 1 сПз

Вязкость нефти 10 сПз

Критические насыщенности воды и нефти 0.25, 1

Начальная водонасыщенность 0.25.

Период разработки 10 лет

Шаг сетки ~23 (м)

Таким образом, имея столь неоднородную модель по геометрическому строению участка, необходимо провести идентификацию интенсивности потоков жидкости между скважинами с помощью «Interwell model», по этим данным определить расположение экрана и трещины.

Постановка задачи для Interwell model. Вернемся к рассматриваемой синтетической модели нефтенасыщенного пласта (рис. 3). В качестве фактических данных будут использоваться объемы добытой жидкости QyK = @нефть + QBoga, закачанной жидкости Q3, среднее пластовое давление Р по скважинам в области их «ареала дренирования» Vi,i = 1...9 (рис. 3). Определим V как цилиндр с центром в месторасположении скважины и с радиусом г = 40 м. Таким образом, поровый объем равен:

^ = пткг2 = п • 0.3 • 10 • 402 « 1.5 • 104 (м3). (9)

Подошвенная вода в настраиваемой модели (2) отсутствует (Апв, £ = 0), вся закачанная вода идет в расчетный пласт (к^ = 1), эффективная сжимаемость постоянна (в £ = 1.25 • 10-5).

Полагаем, что матрица А (3) включает только «сильные» связи, т. е. связи только между добывающей и нагнетательной скважинами, а более слабыми связями между добывающими скважинами мы пренебрегаем. Таким образом, начальный вектор управлений (4) для 9 скважин включает в себя Ыс = 20 связей.

Результаты

Параметрическая идентификация. Проведем параметрическую идентификацию сразу для Ыс = 20 связей. Результаты работы алгоритма параметрической идентификации приведены ниже на рис. 5. Для удобства чтения на рис. 5 б представлены только отфильтрованные значения:

Т- ■ > Т ■

11] — 1тт,

где в качестве порогового значения принималось Тт1П = 10-.

Анализируя рис. 5б, видно, что сильная связь «1-6» может указывать на наличие трещины, а отсутствие связей «1-4», «4-5» показывает, что здесь возможен экран. Однако, есть и некоторые «некорректные» связи: «1-8», «5-8», «4-9», которые не позволяют однозначно говорить о наличии столь протяженного экрана. Значения функционалов (7), (8) представлены в общей табл. 2.

С точки зрения отклонения расчетных значений от фактических (7) найденное значение 1 = и0рг, однако, с точки зрения адекватного отображения действительности результат параметрической идентификации сразу всех связей неудовлетворителен. Необходимо проводить структурную идентификацию.

1пТ,7

а

Ь

Рис. 5. Адаптированные значения Т^ при = 20:

а - показаны все значения; Ь - отфильтрованные Ту > Тт1П

Структурно-параметрическая идентификация. Разработан алгоритм структурной идентификации на основе генетического алгоритма, где в качестве управляющих параметров выбирается вектор тп, определяющий отсутствие или наличие связей между скважинами в матрице А. Элементы вектора тп могут принимать только 2 значения: либо 0, либо 1. Длина

вектора может изменяться в интервале (1,max(Nc)], в нашем случае (1, 20]. Если элемент вектора тп равен 0, то соответствующая связь отсутствует. Формирование матрицы связей А = А'(тп,тп0) включает в себя априорные знания Iapriori = тп0 о наличии/отсутствии связей, в рассматриваемом примере таковыми знаниями будет пренебрежение связями между добывающими скважинами и между нагнетательными скважинами.

Тогда решением задачи структурной идентификации будет такой вектор mn0pt, при котором целевой функционал структурной идентификации минимален

Jstr (u0pt(mn0pt, тп0)) ^ min. Целевой функционал по структурной идентификации может

быть определен формулами (7), (8) или их синтезом.

Применив алгоритм структурной идентификации к нашей задаче, были получены результаты для лучшей с точки зрения функционала модели (рис. 6).

1пТ,7

Рис. 6. Адаптированные значения Т^ при = 9

Статистический подход к структурной идентификации. Также был проведен статистический анализ 10 лучших моделей с точки зрения функционала ]ах х: для каждой связи была получена величина вероятности ее появления в Ыт = 10 лучших моделях:

Tnn-stat =

£ тщ Nm '

где i = 1, ...,Мт обозначает порядковый номер модели.

В качестве критериев отбора моделей используются следующие неравенства:

1. 0 < тщш,) < Рты - ^ой связи нет;

2. Рты < тп3ш,] < Ртах - ^ая связь возможна;

3. Ртах < тщШ1,] <1-}-ая связь есть.

Значения ртхп, ртах выбирались экспертно, Ртп = 0.4, ртах = 0.8. Если выбрать такие критерии отбора моделей, то получается следующее правило отбора (рис. 7): 7 связей из 20 признаются устойчивыми и точно определенными (ртах < тпз1а1^ < 1), 6 связей из 20 определяются как несуществующие (0 < тпз1а1^ < рт{п), про оставшиеся 7 из 20 связей ничего точно сказать нельзя (рт[п < тп51а1^- < ртах). Из рисунка видно, что единственная точно определенная длинная связь - это связь «1-6», где расположена трещина. Область, где располагается «экран», определяется повышенной плотностью отсутствия связей, однако, есть связь «1-8», которая возможна, вероятность ее появления равна 0.6.

Таким образом, за одну итерацию статистического подхода к структурной идентификации были определены 13 из 20 связей, которые отражают действительность (наличие экрана и трещины в модели). Далее, для получения значений этих связей, проведем 2-ю итерацию структурной идентификации, только теперь размерность вектора управляющих параметров тп равна 7, а не 20. Остальные 13 из 20 связей идут в качестве априорной информации тщ.

В результате работы алгоритма структурной идентификации была получена модель (рис. 8).

1пТ,7

Рис. 7. Результат статистической обработки 10 лучших по значению функционала моделей рт1п = ^,ртах = 8

Рис. 8. Адаптированные значения Тц при = 11

Все 11 связей адекватно описывают реальность (рис. 3), что можно считать успешным завершением структурной и параметрической идентификации. Результаты пошаговой структурной идентификации с точки зрения функционалов (7), (8) для каждой из посчитанных моделей (рис. 5, рис. 6, рис. 8) приведены ниже (табл. 2).

Таблица 2

Результаты пошаговой структурной идентификации

№ модели Кол-во связей Ыс В1С За11

1 20 2290 3.24

2 9 2698 3.16

3 11 2392 3.07

Заключение

Для построения фильтрационной модели нефтяного месторождения, позволяющей не только адекватно описывать историю разработки, но и хорошо прогнозировать промысловые показатели, необходима не только параметрическая, но и структурная идентификация. Показано, что для структурной идентификации фильтрационной модели в условиях недостатка априорной информации в целевой функционал, относящийся к истории разработки, требуется вводить дополнительные слагаемые. В настоящей работе предложена методика структурно-параметрической идентификации сеточной модели фильтрации при наличии трещин и экранирующих разломов с использованием только промысловых данных. Методика включа-

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

Литература

1. Мирзаджанзаде А. Х., Хасанов М. М., Бахтизин Р. Н. Моделирование процессов нефтегазодобычи. Нелинейность, неравновесность, неопределенность. М. ; Ижевск : Институт компьютерных исследований, 2004. 368 с.

2. Сергеев В. Л., Наймушин А. Г., Нгуен К. Х. Адаптивная идентификация эволюционных процессов нефтегазодобычи на основе интегрированных систем феноменологических моделей // Изв. Томск. политех. ун-та. Инжиниринг георесурсов. 2014. T. 324. № 1. С. 177-183.

3. Карабутов Н. Н. Структурная идентификация систем: Анализ динамических структур. М. : МГИУ, 2008. 160 с.

4. Каневская Р. Д. Математическое моделирование гидродинамических процессов разработки месторождений углеводородов. М. ; Ижевск : Институт компьютерных исследований, 2002. 140 с.

5. Rafiee M. M. Model Selection and Uniqueness Analysis for Reservoir History Matching : Ph.D. thesis. Freiberg University of Mining and Technology, 2010. 112 р.

6. Косяков В. П., Ширшов Я. В., Мусакаев Э. Н. Вычислительная технология расчета материального баланса на нефтяном месторождении // Нефтепромысловое дело. 2015. № 11. С. 30-35.

7. Льюнг Л. Идентификация систем. Теория для пользователя. М. : Наука, 1991. 432 с.