Математическая модель томографии на давлениях при контроле за разработкой нефтяных месторождений Текст научной статьи по специальности «Физика»

УДК 622.245.5:519.86

МАТЕМАТИЧЕСКАЯ МОДЕЛЬ ТОМОГРАФИИ НА ДАВЛЕНИЯХ ПРИ КОНТРОЛЕ ЗА РАЗРАБОТКОЙ НЕФТЯНЫХ МЕСТОРОЖДЕНИЙ

А.И. КОБРУНОВ

Ухтинский государственный технический университет, г. Ухта [email protected]

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

Ключевые слова: проницаемость, давление, скважины, моделирование распространения, томография, обратные задачи, итерационные процессы, принцип Беллмана, оптимизация, критерий оптимальности

A.I. KOBRUNOV. MATHEMATICAL MODEL OF TOMOGRAPHY AT PRESSURES AT CONTROL OVER THE DEVELOPMENT OF OIL FIELDS

The problem of spatial distribution of effective resistance to fluid movement in permeable formation is considered. The decision method is based on a tomographic processing of time measuring system of reaction approach in wells receivers, on pressure mode change in wells sources. The inverse problem is solved on the basis of the iterative tomographic scheme. At a rare network of data the uniqueness of the decision is provided with implementation of a special design of optimality criterion using a correlation matrix for required parameters.

Key words: permeability, pressure, wells, distribution modeling, tomography, inverse problems, iterative processes, Bellman's principle, optimization, optimality criterion

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

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

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

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

вившегося режима давления в пласте, или обратной ей величине - сопротивлению этому распространению. Этот параметр определяется как время

t прошедшее от подачи давления в интервале пласта в одной скважине, имеющей координаты устья «1 =|x1,у^ , до появления установившегося режима в другой, с координатами «2 = { х2, у2}

отнесенное к расстоянию между I («1,«2) скважинами. Этот параметр меняется по площади, отражая проницаемость пласта в разных точках относительно передачи давления и, следовательно, продуктивности. Здесь влияют меняющиеся коллекторские свойства, глубина залегания, образовавшиеся в процессе движения нефти локальные «заторы» и т.д. В связи с различием проницаемости относительно фаз, установление давления по каждой из них происходит по некоторому монотонному (возрастающему, если подан избыток давления и убывающему, если поданное давление отрицательно -сброс давлений) закону. Рассматривается полностью установившийся режим относительно всех фаз -интегральный стационарный режим.

Формализованная постановка задачи реконструкции скорости распространения давления

Рассмотрим время достижения фронта установившегося давления (далее просто фронта давления) с источником давления в точке «1 в точку

«2. Пусть это время есть т(«1,«2) . Выберем произвольное I > 0 и построим область D(«1,«,t)

занятую установившимся давлением, которая ограничена поверхностью /(«1,«,t) = 0, разделяющую D («1,«, I) и невозмущенную область /(«1,«, I). - это фронт давлений. Если

/(«1,«, I )< 0, то точка « =(X, у) находится внутри D.

В противном случае /(«1,«, I)> 0 и « =(X, у) находится на фронте /(«1,«,I) = 0 , либо впереди него /(«1,«, I)> 0.

Выделим точку «2 на фронте /(«1,«,I) и для произвольной линии L («1,«2) , целиком лежащей в D («1,«, I) , найдем время распространения давления от точки «1 к «2 вдоль Ь («1,«2) по правилу:

г(«1’«2)= I Й7Й= I к(х>у)^1

Ь(«1 ,«2 ) У у) Ь(«1 ,«2 )

(1)

Здесь V ( X, у) - скорость, а к (X, у) - обратная ей «медленность» передачи давления в точке (X, у) . Считается, что эта величина неоднородна, но изотропна. Исследуемые функции ограничены и интегрируемы в силу конечности всех рассматриваемых областей. Впервые давление

придет в точку «2, двигаясь по траектории, обеспечивающей минимальное из всех возможных времен и, следовательно, это первое время I есть минимальное из всех возможных по разным траекториям:

= min

Ь(«1,«2 )

аг

Щ,«2 )

V ( X, у )

(2)

Соотношение (2) - это аналог интеграла Ферма, известного в геометрической оптике [1] и теории распространения сейсмических и электромагнитных волн, распространенного на задачу о движении давлений. Считая траектории I(«1,«2) заданными, определим оператор:

АИ«)] = I««,), (3)

отображающий распределение скорости передачи давления, сгенерированного в скважине, с координатами « во времена первых вступлений давления в скважину, находящуюся в «, . Строго говоря, этот оператор следовало бы записывать с учетом значения индексов /, ] . Однако имея в виду соотношение (2), мы опускаем их, считая вполне достаточным указание на их значения в правой части (3). Рассмотрим обратную для (3) задачу реконструкции пространственной скорости распространения давлений в пласте, по известным временам достижения установившегося режима по всем парам скважин. Число таких пар - N . Следовательно, размерность входных параметров также N . Для ее решения воспользуемся приемом линеаризации.

Пусть задано начальное V0 («) приближение

к полю скоростей V («) . И пусть поправка

5(X,у) = V(X,у)-V0 (X,у) .

(4)

Тогда:

А [V («)] = А \у0 («)] + А'[У0 («)>(«)+г (5)

. (5)

Здесь г (5) малая высшего порядка относи-

тельно второго члена в

(5), Л'[Г° (Ц)]

м

Ут (Ц) = д

производ-

(11)

ная Фреше [2] для (3). Легко вычислить, следуя традиционной лучевой томографической технологии, см. например [3] пренебрегая малыми высшего порядка в сравнении с 5(«) :

А'|>0 («)]5(«) = А [V" («) + 5(«)] —А [V0 («)] + г (5):

ё! г ё!

Г ш _ Г „¿.)У°(*-У) + 5(И) «і,)У’ (х,У)

+ Г

(«),

/(ЙА )

[у0 (х,у)_{у0 (х,У) + 5(1)}]С §(Ц)сЧ

[У0 (х,у) + <5(Ц)]-У0 (х,у)

У2 (и) (6)

Тогда,пренебрегая членом г(5):

ёг г 5 («) аг

Л \У (и)]= I

I

Уравнение (3) перепишется: 5(Ц ) С/

(

( ¿)(у 0 )2 (и )

(7)

(8)

где:

д( Ц Ц, ) = і0 (Ц Ц,)_ і (Ц,, Ц,);

С1

(9)

У0 (х,У)■

Алгоритм обращения

Рассмотрим задачу нахождения решения 5( « ) уравнения (8) по заданным правым частям

Д(«i, « ,) , образующим в своей совокупности вектор А = {Дп, п = 1,...N} , N - число доступных

для анализа пар и {/,у, / Ф,} . В соответствии с технологией проекционных методов [4] будем искать 5(« ) в виде представления:

м

5(И ) = Е Чт Ут (Ц )

(10)

т=1

где /т ( « ) - заданная базисная система функций. Например, это постоянные функции в непере-секающихся подобластях Sm , в совокупности своей покрывающих всю изучаемую площадь. Подставив (10) в (8) и учитывая выражение для А полу-

т=1 4 (V0 ) ( « )

п = 1 * N.

Здесь гп = I («г,«,) , где («г,« ,) образуют

пару с номером п .

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

а = Кп};

// ( « ) аг

а..

(12)

(У0 )2 (и)

Следует найти решение системы линейных уравнений:

ая = А;

Г !, (13)

Я = {Чт , т = 1 *М}

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

Обращение матрицы (13). Рассматривая вектора я и А как элементы евклидовых про-

Г)М T)N

странств К и К соответственно, заменим задачу (13) на:

||ая_Д||2 +а F\я_q0]

(14)

Здесь F невырожденная М хМ матрица, например единичная. Она определяет вид критерия

оптимальности

г[я_я°]

характеризующего

чим:

меру уклонения искомого вектора Я и нулевого

приближения к нему Я0 . Например, представляет собой обратную к корреляционной матрице для

[я — Я0 ] . Вопросы выбора этой матрицы относятся к методическим вопросам нахождения оптимальных решений [5]. О - числовой параметр, параметр регуляризации задачи, определяющий точность решения уравнения в (13). Значение оптимального параметра регуляризации подбирается по принципу невязки [5]. Варьируя квадратичный функционал (14), получим для ее решения [3,5] уравнение (уравнение Эйлера для вариационной задачи) в виде:

а* (ая — А) + аР*Р [я — Я0 ] = 0. (15)

Звездочка над операторами ( матрицами), как обычно, означает переход к сопряженным (в данном случае транспонированным) операторам. Задача (15), в

п

2

2

отличие от (13), состоит в обращении уравнения с квадратной и невырожденной матрицей:

[а*а+aF*F ] я = а*А+aF*Fя0 (16)

Для этой цели могут быть применены стандартные алгоритмы обращения:

Я(а) = [а*а + aF*F] [а* А+aF*Fя0] . (17)

Расчет траекторий 1п . Вычисление элементов атп матрицы а по формуле (13) требует знания траекторий 1п , соединяющих по принципу минимального времени (2) пару точек !п = I («.,«.) . Основой их нахождения служит анализ поведения линий фронтов /(«1, «, I ) = 0 во времени и использование принципа оптимальности Беллмана [6] для нахождения траектории, реализующей минимальное время движения в паре («г, «1 ) . Зафиксируем точки («г, « 1 ) . Время I («г, « 1 ) распространения давления из «г в «. определено как решение уравнения /(«г,«,, I ) = 0 относительно

I. Принцип Беллмана [6], примененный к задаче построения оптимальной траектории (луча)

I(«1,« 1) , состоит в том, что если п - произвольная точка, лежащая на г (« «1), то г (п,« 1) есть

оптимальная траектория (луч) из точки п в «..

Иными словами, каждый завершающий фрагмент оптимальной траектории - суть самостоятельная оптимальная траектория, соединяющая соответствующие точки. Этот принцип может быть сформулирован и иначе. Оптимальность траектории от промежуточной точки не зависит от того, как луч попал в эту точку. Распространение луча - есть система без памяти. Алгоритм Беллмана численного дискретного построения траектории состоит в следующем:

1. Интервал отрезка, соединяющего точки «1 и «1, разобьём на Nx частей длиной ДЬХ . Номера интервалов гх = 1,2,...Nx . В каждом узле проведем ортогонально отрезку прямую, разбив ее на интервалы ДЬТ . Число этих интервалов NY.

Номера интервалов iY = 1,2,...NY . Таким образом, организована локальная координатная сетка с точка-

(^, у*). (-V у0 ) = «г; (Х, у^) =«1.

ми

<■" (xn,) = t [(xn,-1,У,,), £,]

ния давления из точек (XN^ ) в £

распростране-Какой-

то их этих

N

отрезков, соединяющих

( XNX -1> угг ), и £j , принадлежит I (£,, £, ) .

3. Далее выделим линию, отвечающую точке

X

N,-2 •

Для каждого х„ -2,у w найдем то у.

из линии с XN -., которое обеспечивает:

min

xNx -1> У'г

ІПу,г ^ [( XNx-2, Уіу ) , ( XNx-1, Угг )] + * ( XNx-1 )}

( XNx-2 )

= Г І X:

(18)

Точка минимума определит новую точку из г( «¡,«1) . Минимизация идет по у ^ из линии с

X

N,-1 •

Одновременно последовательно наращиваются фрагменты вариантов траектории, общее число которых всегда остается Ыт . Получим в результате Ыт значений г'т (хм^_, ) .

4. Продолжим процесс аналогичным обра-

зом:

{' [(X,X ■ Угг ) ■(X,.!• Угг )] + (\ )}= ''' (X'x ) ■

х'Угг v ’

(19)

В конечном итоге будет достигнута точка

и вычислено N значений:

X

2. Выделим линию, отвечающую точке вектор t = {, = ,£ £ ) ' = . ^n} '

( Уо ) = £ '

,Y ( X0 )= min {t \ (£' ), ( X,+1> Уг )] + ,Y ( X1 )}■

X'x ,

(20)

5. Простым выбором минимума из Nr значений будет определена истинная траектория луча

1 (£'’£j) ■

Алгоритм уточнения пространственной модели эффективной скорости распространения установившегося режима давлений - томографических изображений состоит в следующем:

Исходные данные:

Координаты £г источника давлений; £, -точки измерения режима давления; t ( £г, £ , ) - регистрируемое время распространения установившегося давления сгенерированного в скважине £г до скважины £ ,, представляющие собой

= 1

N, -1

. Подсчитаем NY времен

Модель нулевого приближения V0 (X, У) ;

1 = 0

Шаг 1. На основе принципа Гюйгенса [1] для каждого «г строится фронт давлений // («1, «, I) ;

Шаг 2. По описанной выше схеме находится семейство траекторий 1'(«г • «1) из точки «г в « 1

для всех г ф 1, г, 1 = 1 * N.

Шаг 3. Вычисляется:

А'={Д,(«, «1 )=< «1)—< («, «1)};

<,()=-(!«1) ■

Шаг 4. Выбираются аппроксимационные функции: /т ( « );

Шаг 5. Рассчитывается матрица:

а1 =

{атп } ;

/т ( « )

I

(V') («)

Шаг 6. Выбрав F , рассчитываем:

Я'+1 (а) = [а*а + aF*F] [а* А + aF*Fя'

М

Шаг 7. 5'+1 (« ) = ]Г / (« ) .

Шаг 8. V1+1 ( X, у) = V1 ( X, у) + 51+1 ( X, у) .

Шаг 9. Контроль точности подбора, числа сделанных итераций, динамики сходимости. В случае необходимости ' = ' +1; Переход к шагу 1. В противном случае КОНЕЦ.

Заключение

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

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

Литература

1. Кравцов ЮА., Орлов Ю.И. Геометрическая оптика неоднородных сред. М.: Наука, 1980. 304 с.

2. Иоффе Ф.Д., Тихомиров В.М. Теория экстремальных задач. М.: Наука, 1974. 479 с.

3. Кобрунов А.И. Математические основы теории интерпретации геофизических данных: учебное пособие. М.: ЦентрЛитНефтеГаз, 2008. 286 с.

4. Марчук Г.И. Методы вычислительной математики. М.: Наука, 1980. 534 с.

5. Некорректные задачи математической физики и анализа / М.М. Лаврентьев и др. М.: Наука, 1980. 286 с.

6. Болтянский В.Г. Математические методы оптимального управления. М.: Наука, 1969. С. 408.

Статья поступила в редакцию 11.04.2012.