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

Виталий Кунцев

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

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

Информатика

Диссертации

Вуз: Ухтинский государственный технический университет (УГТУ)

ID: 5f4cf159cd3d3e0001ce9658

UUID: 08c82440-cdb6-0138-26db-0242ac180006

Язык: Русский

Опубликовано: больше 3 лет назад

Просмотры: 9

10.18

Виталий Кунцев

Комментировать 0

Рецензировать 0

Скачать - 12,3 МБ

Поделиться работой

МИНОБРНАУКИ РОССИИ Федеральное государственное бюджетное образовательное учреждение высшего образования «Ухтинский государственный технический университет» (УГТУ) На правах рукописи Кунцев Виталий Евгеньевич РАЗВИТИЕ МЕТОДА ГИДРОДИНАМИЧЕСКОЙ ТОМОГРАФИИ НА ОСНОВЕ СИНТЕЗА ДАННЫХ ИЗ МАТЕМАТИЧЕСКОЙ МОДЕЛИ ЭКСПЛУАТАЦИИ НЕФТЕГАЗОВОГО МЕСТОРОЖДЕНИЯ 05.13.18 – Математическое моделирование, численные методы и комплексы программ Диссертация на соискание ученой степени кандидата технических наук Научный руководитель: доктор физико-математических наук, профессор, А. И. Кобрунов Ухта – 2018

ОГЛАВЛЕНИЕ Введение ....................................................................................................................... 4 Глава 1. Обзор методов реконструкции фильтрационного сопротивления проницаемых пластов ............................................................................................... 11 1.1. Гидродинамическое исследование проницаемых пластов....................... 11 1.2. Фильтрационно-емкостные свойства проницаемых пластов .................. 18 1.3. Томографическая система гидропрослушивания скважин ...................... 19 1.4. Математическое моделирование эксплуатации скважин нефтегазового месторождений .......................................................................................................... 21 Глава 2. Пассивная гидродинамическая томография ............................................ 26 2.1. Математическая модель динамики эксплуатации месторождения ......... 26 2.2. Итерационный алгоритм подбора оптимального вектора параметров ... 35 2.3. Заключение .................................................................................................... 43 Глава 3. Методика и технология решения задачи гидродинамической томографии проницаемого пласта ........................................................................... 44 3.1. Организация и обработка результатов гидродинамического прослушивания скважин в томографическую систему данных ........................... 44 3.2. Алгоритм расчета интервальных времен движения особых точек динамики восстановления давления........................................................................ 46 3.3. Вычислительная схема метода гидродинамической томографии ........... 55 3.4. Апробация вычислительной схемы гидродинамической томографии ... 59 3.5. Заключение .................................................................................................... 65 Глава 4. методы оценки результатов работы гидродинамической томографии 67 4.1. Построение атрибута нечеткой меры по результатам гидродинамической томографии ............................................................................. 67 2

4.2. Построение нечеткой модели распределения фильтрационного сопротивления проницаемого пласта ...................................................................... 77 4.3. Заключение .................................................................................................... 85 Глава 5. Программный комплекс моделирования пассивной гидродинамической томографии проницаемого пласта ........................................................................... 87 5.1. Функциональные возможности ................................................................... 87 5.2. Блок-схемы алгоритмов пассивной гидродинамической томографии ... 89 5.3. Пользовательский интерфейс программного комплекса ......................... 91 5.4. Заключение .................................................................................................... 98 Заключение................................................................................................................. 99 Список литературы.................................................................................................. 101 Приложение А.......................................................................................................... 114 Приложение Б .......................................................................................................... 121 3

ВВЕДЕНИЕ Актуальность темы исследования. На сегодняшний день существует проблема расширения области применимости методов изучения объектов с использованием томографических принципов, например, в виде получения послойных изображений внутренней структуры объекта, поскольку применение методов томографии требует адаптации и доработки конкретных модификаций методов под задачи, на которые они расширяются. В этой связи большим потенциалом обладают задачи нефтегазопромысловой отрасли, где томографическая обработка позволит более детально изучить исследуемый объект, в частности с помощью гидродинамической томографии изучить фильтрационно-емкостные свойства месторождения. При решении подобных задач томографическим методом, основная проблема томографии – это проблема «размерности» задачи, т. е. необходимо большое количество исходных данных, для получения которых необходим определенный комплекс технологических работ. В условиях высокой экономической и временной себестоимости процесса получения этих данных, метод гидродинамической томографии исследования фильтрационно- емкостных свойств оказывается слишком дорогостоящим, поэтому была предложена идея научиться синтезировать необходимые исходные данные, используя методы математического моделирования. Задача построения достоверной фильтрационно-емкостной модели нефтегазового месторождения является актуальной, так как часто возникает проблема ухудшения фильтрационных свойств проницаемого пласта, особенно на поздних стадиях эксплуатации месторождений, что может выражаться в понижении уровня добычи скважин. В связи с этим проблема борьбы с эффектами и результатами ухудшения фильтрационных свойств пласта является очень актуальной. Основным компонентом решения данной проблемы является точное определение местоположения в пространстве участков аномального фильтрационного сопротивления. 4

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

принципов. Математическое синтезирование этих данных возможно по результатам изучения истории динамики добычи скважин на временных интервалах штатной эксплуатации нефтегазового месторождения, после которых стало существенно выражаться понижение фильтрационного сопротивления в виде уменьшения объемов добываемых углеводородов. Для решения задачи необходимо разработать модель прогнозирования динамики добычи месторождения по истории эксплуатации, отражающей изменение пространственных характеристик фильтрационного сопротивления проницаемого пласта, с помощью которой планируется проведение гидродинамического прослушивания скважин и преобразование полученных данных в томографическую систему. Степень разработанности темы исследования. Задачам исследования пространственного распределения фильтрационных параметров проницаемых пластов, характеризующих пропускную способность продуктивных пластов посвятили своим работ многие ученые: Щелкачев В. Н., Чарный И. А., Умрихин И. Д., Бузинов С. Н., Лейбензон Л. С., Борисов Ю. П., Баренблатт Г. И., Сургучев М. Л., Басниев К. С., Каневская Р. Д. Карнаухов В. Л., Пьянкова Е. М., Кременецкий М. И., Ипатов А. И., Мирзанджанзаде А. Х., Gringarten А. Joshi S., Raghavan R и др. Построением моделей связности скважин занимались: Краснов В. А., Иванов В. А., Хасанов М. М., Valko P. P., Doublet L. E., Yousef, Jong S. Kim, Larry W. Lake, Thomas F. Edgar. В научных работах Кобрунова А. И., Дорогобед А. Н. Куделина С. Г. описана общая технология проведения технологических мероприятий для реализации гидродинамической томографии для действующих скважин нефтегазового месторождения. Тем не менее, существует ряд проблем, требующих решения. В следствие этого необходимо развитие новых методов решения этих проблем. Основным отличием работы является непрерывный мониторинг динамики изменения эксплуатационных характеристик с целью прогноза текущего состояния проницаемого пласта и появления новых зон с аномальными фильтрационными характеристиками. 6

Цель диссертационного исследования. Разработка математических методов моделирования, основанных на синтезе данных из пассивной и активной форм гидродинамической томографии для обнаружения и локализации пространственных зон нарушения фильтрационных характеристик проницаемого пласта, построение численных алгоритмов, соответствующих этим методам, и их компьютерная реализация. Основные задачи диссертационного исследования. 1. Построить математическую модель динамики эксплуатации месторождения с целью синтезирования томографических данных (пассивная гидродинамическая томография). 2. Разработать вычислительную схему, прогнозирования пространственного распределения фильтрационных параметров на основе решения обратной задачи, адаптированную к синтезированным данным пассивной гидродинамической томографии. 3. Разработать метод оценки результатов прогноза пространственного распределения фильтрационных параметров. 4. Разработать комплекс прикладных программ для реализации математической модели и вычислительных схем пассивной гидродинамической томографии с целью прогнозирования пространственного распределения фильтрационных параметров. Объект исследования: фильтрационно-емкостные модели нефтегазовых месторождений и история эксплуатации скважин. Предмет фильтрационного исследования: сопротивления пространственное проницаемого пласта распределение нефтегазовых месторождений и правила оценки результатов поиска аномальных зон фильтрационного сопротивления. Методы исследований. В диссертационном исследовании для решения поставленных задач используются методы математического моделирования, решения обратных задач, элементы функционального анализа, динамического программирования и теории нечетких множеств. Для проектирования и 7

реализации программного комплекса использовались методы объектноориентированного программирования. Научная новизна диссертационного исследования. 1. Выполнено математическое моделирование динамики работы месторождения, позволяющее синтезировать данные для использования схем реконструкции фильтрационных параметров. 2. данных, Разработана вычислительная схема обработки томографических позволяющая фильтрационных прогнозировать параметров на пространственное основе решения распределение обратной задачи гидродинамической томографии. 3. Разработан математический метод, основанный на использовании интервальных оценок для исходных данных, позволяющий получить оценку результатов прогноза пространственного распределения фильтрационных параметров. 4. Разработан комплекс прикладных программ для реализации математической задачи гидродинамической и томографии вычислительных с схем проведением пассивной вычислительных экспериментов. Теоретическая и практическая значимость. Разработаны методы и алгоритмы пассивной гидродинамической томографии для прогнозирования пространственного распределения фильтрационного сопротивления для локализации зон аномального фильтрационного сопротивления с целью дальнейшего проведения необходимых геолого-технических мероприятий (ГТМ) по ликвидации застойных участков в межскважинном пространстве проницаемых пластов нефтегазовых месторождений. Метод отличается тем, что с целью снижения экономических и временных затрат на проведение работ по проведению гидродинамического прослушивания скважин, исходные данные для томографической обработки рассчитываются по математической модели работы месторождения исходя из динамики его эксплуатации. 8

Результаты, полученные при выполнении диссертационного исследования, применяются в учебном процессе в Ухтинском государственном техническом университете на кафедре геофизические методы, геоинформационные технологии и системы. Положения, выносимые на защиту. Математической метод моделирования и вычислительная схема 1. пассивной гидродинамической томографии для синтезирования данных с целью включения их в активную форму томографии. Математический 2. метод и алгоритм прогнозирования пространственного распределения фильтрационных параметров на основе решения обратной задачи томографии с использованием синтезированных данных. Математический 3. метод оценки результатов прогноза пространственного распределения фильтрационных параметров с учетом погрешности исходных данных. Комплекс прикладных программ для реализации математической 4. задачи и вычислительных схем пассивной гидродинамической томографии с целью прогнозирования пространственного распределения фильтрационных параметров и анализа результатов. Степень достоверности и апробация результатов. Апробация выполнялась на фактическом материале по одному из месторождений, относящихся к Тимано-Печорской нефтегазоносной провинции. Результаты апробации показали возможность получения необходимых интервальных времен на основе разработанных алгоритмов и построенной модели эксплуатации скважин нефтегазового месторождения. Диссертационная работа прошла внедрение в научно-исследовательском центре Коми регионального отделения Российской академии естественных наук «Институт геотехнологий». Основные положения и результаты диссертации были представлены и обсуждены на международной молодежной 9 научной конференции

«Севергеоэкотех» (Ухта, 2015), VI Международной студенческой электронной научной конференции «Студенческий научный форум 2015», международном семинаре «Рассохинские чтения» (Ухта, 2015-2018), международном семинаре им. Д. Г. Успенского (Пермь-2015, Воронеж-2016, Москва-2017, Казань-2018), III Школе-конференции «Гординские чтения» (Москва, 2015), республиканском молодежном инновационномконвенте «Молодежь – будущему Республики Коми» (Ухта-2015), 37-й Всероссийской научно-практической конференции с международным участием «Геология и полезные ископаемые Западного Урала» (Пермь, 2017), XVIII Международном семинаре «Физико-математическое моделирование систем» (2017, Воронеж), а также на научных семинарах кафедры информатики, компьютерных технологий и инженерной графики и кафедры геофизических методов, геоинформационных технологий и систем Ухтинского государственного технического университета. По материалам диссертации опубликовано девятнадцать работ, из которых пять – в журналах из Перечня российских рецензируемых научных журналов, четырнадцать статей – в материалах конференций. На программный государственной «Прогнозирование комплекс регистрации работы получены программы скважин для два свидетельства ЭВМ: нефтегазового о № 2017619300 месторождения», № 2017662707 «Пассивная гидродинамическая томография проницаемого пласта» и № 201817703 «Построение интервальных оценок гидродинамической томографии проницаемого пласта». Структура и объем диссертационной работы. Диссертационная работа состоит из введения, пяти глав, в которых отражены основные результаты исследования, заключения, списка использованных источников литературы и приложения. Общий объем диссертации составляет 141 страницы с 62 рисунками. Список литературы включает 111 наименований. 10

ГЛАВА 1. ОБЗОР МЕТОДОВ РЕКОНСТРУКЦИИ ФИЛЬТРАЦИОННОГО СОПРОТИВЛЕНИЯ ПРОНИЦАЕМЫХ ПЛАСТОВ 1.1. Гидродинамическое исследование проницаемых пластов Одним из способов, который может помочь улучшить эффективность процесса разработки месторождения [60], является построение достоверной фильтрационно-емкостной модели нефтегазового месторождения. Но построение адекватной фильтрационно-емкостной модели продуктивных коллекторов эксплуатируемых месторождений труднореализуемо из-за неизвестности фильтрационно-емкостных характеристик в межскважинном пространстве месторождения, кроме этого в процессе разработки ввиду движения флюидов внутри проницаемого пласта происходят различные химические процессы, которые могут существенно повлиять на ранее выявленные фильтрационно-емкостные свойства пласта. Происходят различные эффекты формирования нефтяных тромбов, участков повышенной вязкости, возникают эффекты асфальтизации и поровое пространство забивается асфальтитами, которые не припускают потоки флюидов. В связи с этим фильтрационная модель меняется, ранее выявленные фильтрационные потоки существенно изменяются, и сетка скважин месторождения больше не может считаться оптимальной под текущие фильтрационно-емкостные свойства пласта. Гидродинамические факторы основаны на гидромеханическом засорении фильтрующих каналов проницаемых пластов различными механическими примесями, которые содержатся в технической воде, закачиваемый в пласт через нагнетательные скважины месторождения, и углеводородными химическими соединениями, образующимися в процессе движения флюидов внутри пласта. К ним могут относиться оксид железа, соли и карбонаты, сульфаты и хлориды, а также частицы песка и глины. К примеру, на глубокое загрязнение порового пространства проницаемого пласта содействует высокое 11

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

зоне, а также в удаленных участках пласта. Сегодня к часто используемым методам ГДИС относят метод кривой восстановления давления (КВД) и метод гидродинамического прослушивания скважин [3]. Гидропрослушивание [63, 69] скважин нефтегазового месторождения выполняется для анализа фильтрационных характеристик проницаемого пласта (пористость, проницаемость, пьезопроводность и др.), изменений в тектонике и т. п. Основной идеей методов гидродинамического прослушивания скважин является мониторинг за изменениями давлений в добывающих скважинах (скважина-приемник), обусловленным изменением режима работы окружающих возмущающих скважинах (скважина-источник). Регистрируя время начала остановки работы или изменения режима отбора флюидов в возмущающей скважине-источнике и время начала изменения режима давления в реагирующей скважине-приемнике, по интервалу времени в течении которого пробегает волна давления внутри пласта между парой рассматриваемых скважин можно проводить анализ полученной информации и делать заключение о фильтрационных свойствах проницаемого пласта в пространстве между скважинами изучаемого месторождения. Дойдя до простаивавшей (или работавшей на постоянном режиме давления) скважины, эта волна повышает давление на участке вокруг скважины. Возможен случай, что при гидродинамическом прослушивании в скважине не происходит реагирование на изменение режима работы в соседней скважине. Отсутствие реакции на повышение давления указывает отсутствие или очень плохую гидродинамическую связь между парой рассматриваемых скважинам из-за образовавшегося непроницаемого экрана (тектоническое нарушение). Таким образом, гидродинамическое прослушивание скважин помогает установить фильтрационные особенности структуры проницаемого пласта, которые не всегда есть возможность геофизического несколько изучения различных отличающихся установить нефтегазового способов технологическими в процессе геологического месторождения. проведения деталями [3]. 13 и Существуют гидропрослушивания, Исследование методом

интерференции возмущением, называется т. скважине-источнике, е. гидропрослушиванием происходит которое длительное вызывает с изменение изменение режима однократным дебита в работы в наблюдаемой реагирующей скважине. В результате интерпретации результатов такого гидропрослушивания можно определить гидропроводность пласта. Импульсное исследование называется гидропрослушивание с многократным возмущением, которое проводится путем создания в возмущающей скважине кратковременных импульсов давления (с менее заметными изменениями давления в реагирующей скважине). Могут быть определены такие параметры пласта, как пьезопроводность, гидропроводность и упругоемкость пласта. Использование результатов гидродинамического прослушивания заключается в проведении анализа динамики изменения давления в локальной окрестности исследуемой скважины на изменение работы в возмущающей скважине и последующую интерполяцию результатов такого исследования по всем скважинам межскважинное пространство месторождения. Также В процессе интерполяции результатов гидропрослушивания дополнительно может быть использована имеющая геолого-геофизическая информация по исследуемому объекту. По определенной дополнительно текущей промыслово-технологической информации о работе каждой скважины строят карты текущей насыщенности агентом вытеснения. С учетом построенных карт математически моделируют процесс фильтрации в слоисто-неоднородном пласте. В качестве дополнительной информации определяют текущее содержание в продукции каждой скважины исследуемого участка агента вытеснения [4]. Для уточнения результатов гидропрослушивания скважин выполняется анализ матрицы корреляций по данным дебитов скважин в виде объемов извлекаемых углеводородов с данными по объемам нагнетаемой воды [5]. Сущность изобретения: выбирают участок нефтяной залежи с добывающими и нагнетательными скважинами. Учитывают для каждой добывающей скважины данные по объемам сбора нефти, жидкости и воды, времени эксплуатации скважины. Согласно изобретения для каждой нагнетательной скважины производят сбор данных по 14

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

На реализацию метода гидродинамического прослушиванию накладывается несколько условий, которые могут существенно затруднить его применение на месторождениях в режиме штатной эксплуатации. Также метод гидропрослушивания и его аналоги, не позволяют выявить локальные пространственные нарушения области, поскольку в проницаемости используемых пласта начальных в межскважинной данных отсутствует информация о значениях фильтрационного сопротивления во внутренних участках пласта между скважинами в пределах нефтегазового месторождения, необходимая и достаточная для однозначного толкования характера пространственного распределения коэффициента проницаемости. Автор работы [7] описывает следующий метод, отличающийся от обычных подходов в проведении гидродинамического прослушивания скважин. В описываемом подходе предлагается оценивать ФЕС межскважинного пространства локального участка нефтеносного коллектора (ЛУНК) путем анализа этих локальных зон с определением параметров порождающих уравнений. Для каждой рассматриваемой скважины указывается своя зона модели локального участка нефтеносного коллектора. Для фильтрационноемкостных параметров центральной зоны участка определяются уравнения [7]. Определение параметров полученного уравнения осуществляется с помощью традиционных методов, в частности, с применением метода наименьших квадратов (МНК) и модели, приводимой его к регрессионному виду. Для идентификации фильтрационного-емкостных свойств внутри межскважинного пространства оперируют некоторыми осредненными значениями переменных, которые недоступны для измерения. Следующим способом нахождения фильтрационных характеристик пласта по результатам наблюдений за системой: закачка – дебит по всей совокупности скважин, служит принцип подбора параметров фильтрационной модели на основании моделирования фильтрационных потоков [8], [9]. Описанные вычислительные технологии, реализующие наиболее общие модели фильтрации, требуют использования многопроцессорных вычислительных 16

систем с распараллеливанием вычислений [10].Однако, эта и подобные им технологии, позволяя в достаточно общем виде эффективно решать прямую задачу – моделировать фильтрационные потоки в самых общих предположениях о параметрах среды и движущейся смеси имеют малые возможности для подбора фильтрационной модели с целью нахождения ее параметров по наблюдаемым данным (дебит–закачка). Эта, последняя задача по своему духу относится к обратным геофизическим задачам и требует специализированных подходов для своего решения. Они должны ориентироваться на построения изображения локальных неоднородностей относительно аномальных фильтрационных характеристик пласта, обеспечивающих их обнаружение и пространственную локализацию, возможно и в ущерб точности оценки величины фильтрационного сопротивления. В работе [11] показано применение полуаналитической модели для поиска и уточнения фильтрационно-емкостных свойств проницаемого пласта пласта на основе штатных промысловых данных по истории динамики дебитов добывающих скважин. В основе описания распределения песчаных тел по размерам лежит вероятностный подход. Характеристики размеров песчаных тел описываются с помощью двухпараметрического распределения. На основе полученного распределения (песчаных тел) моделируется работа скважины, вскрывающей расчлененный пласт. Описание работы данной скважины сводится к моделированию многопластовой системы с использованием сопряжения решений уравнения пьезопроводности (с различными граничными условиями) на граничных условиях на скважине, причем геометрические размеры пластов вычисляются из параметров распределения песчаных тел. Скорость расчета обеспечивается относительной простотой построенной полуаналитической модели (например, с расчетом на гидродинамическом симуляторе) и позволяет ее использовать для решения некоторого класса обратных задач, в частности, для определения характеристик пласта по истории динамики дебита скважин, что делает возможным ее применение для уточнения 17

фильтрационно-емкостных характеристик проницаемых пластов нефтегазовых месторождений. 1.2. Фильтрационно-емкостные свойства проницаемых пластов Фильтрационно-емкостные свойства (ФЕС) определяют способность проницамых пластов вмещать и накапливать углеводороды (пористость) и фильтровать (проницаемость) флюиды при изменении давления. Под пористостью [13] коллектора определяется наличие в нем мелких пор, трещин и других свободных полостей, вмещающих флюиды, в виде нефти, газа и воды. Коэффициент пористости определяет ёмкостные свойства продуктивного коллектора. Проницаемость - это свойство породы пропускать жидкость или газ при перепаде давления. Проницаемость зависит от размеров и формы поровых каналов. Единицей измерения проницаемости является Дарси [13]. Фильтрационно-емкостные свойства количественно могут быть определены через коэффициент пъезопроводности [15], пространственное распределение которого детально описывает пространственное распределение фильтрационного сопротивления проницаемого пласта. Харченко Б. С. в 1957 году впервые предложил находить коэффициент пьезопроводности проницаемого пласта по итогам анализа кривой восстановления давления. Позже в 1984 году данная задача была развита Г. И. Баренблантом [16, 18], Ю. П. Борисовым [17-22] и В. Н. Щелкачевым [23, 24, 59], И. А. Чарный [58]. В работах представленных авторов была развита методика анализа данных, извлекаемых из добывающей скважины, которая находится под эффектом работы в нагнетательной скважине. Результатом такой методики является только одно значение коэффициента пьезопроводности, который в дальнейшем рассматривается как фильтрационно-емкостная характеристика всего участка межскважинного пространства между парой рассматриваемых скважин. В настоящее время фильтрационно-емкостные свойства проницаемого пласта могут быть определены путем рассматорения 18

неустановившихся процессов, происходящих в пласте при остановке и запуске отдельных скважин [29]. 1.3. Томографическая система гидропрослушивания скважин Возможность изучения пространственного изображения неоднородного (содержащего локальные аномальные минимумы и максимумы) распределения проницаемости пласта в пределах разбуренного участка месторождения в процессе его эксплуатации существует и основана на специализированной методике изучения интервальных времен распространения установившегося режима давления между парами скважин и современных информационных технологиях томографического типа [25, 64]. Задача томографической обработки проницаемого пласта состоит в том, чтобы представить достоверную информацию о пространственном распределении переменной эффективного фильтрационного сопротивления, имеющей характер пропускной способности флюидов внутри пласта под воздействием стационарного режима давления по площади, с целью информационной поддержки технологических мероприятий по увеличению объемов добычи на месторождении путем обнаружения и ликвидации найденных зон проницаемости пониженной пласта, проводимости выводящего или полного фрагменты купирования месторождения из эксплуатации. Способ гидродинамической томографии проницаемого пласта [26, 27] направлен на нахождение томографического изображения проницаемости пласта, основан на специализированном формировании томографического набора данных для распространения флюида в пласте и последующего построения изображения с использованием томографической модификации методов интегральной геометрии, адаптированной к кинематическим уравнениям движения установившегося потока флюидов в неоднородной среде с использованием оптимизационных принципов. 19

Основными методов элементами интегральной активной геометрии, томографической адаптированной к модификации кинематическим уравнениям движения установившегося потока флюидов в неоднородной среде с использованием оптимизационных принципов, служат: 1. построение кинематической модели движения установившегося потока флюидов в виде модели конечных элементов для функции скорости движения потока [28]; 2. моделирование интервальных времен движения установившегося режима движения флюидов по всем парам скважин источник - приемник месторождения на основе оптимизационных принципов [29]; 3. расчет невязки между реально измеренными (экспериментальными интервальными временами) и смоделированными интервальными временами; 4. решение обратной задачи интегральной геометрии [31] для рассчитанной невязки и нахождение поправки к изначальной модели; 5. контроль точности полученной модели пространственного распределения фильтрационного сопротивления проницаемого пласта. Необходимые первоначальные данные по интервальным временам пробега волны давления между парами скважин для осуществления томографического исследования и обработки информации берутся из прямого эксперимента на остановленном месторождении – возникает активная модификация гидродинамической томографии. Это достаточно затратная и технологически трудоемкая схема, которая включает в себя следующие этапы: 1. создание стационарного аномального (положительного либо отрицательного в зависимости от условий эксплуатации месторождения) давления в скважине-источнике аномального давления в интервале продуктивного пласта; 2. регистрацию времен распространения установившегося режима давления в том же участке пласте от скважины-источника до всех скважин, которые служат приемниками, распределенных по месторождения образующих верную систему наблюдений; 20 сети на площади

3. последовательное повторение верной регистрации интервального времени распространения установившегося давления с использованием в качестве скважин источника аномального давления всех скважин, участвующих в веерной регистрации; 4. обработкой полученных двухиндексных данных интервального времени по томографической модификации методов интегральной геометрии, адаптированной к кинематическим уравнениям движения установившегося потока флюидов в неоднородной среде с использованием оптимизационных принципов. 1.4. Математическое моделирование эксплуатации скважин нефтегазового месторождений Математическое описание процессов разработки нефтяных залежей [5153] имеет своей целью предсказание локальных или интегральных характеристик пластовой системы при различных условиях воздействия на нее, нахождение оптимальных режимов эксплуатации нефтяных месторождений. «Для математического моделирования процессов нефтедобычи используются различные подходы, где основные результаты получены путем гидродинамического описания и анализа динамики пластовой системы с помощью дифференциальных уравнений, решения возникающих краевых задач. Такое описание опирается, обычно, на информацию, получаемую при выполнении поисковых и разведочных работ, бурения и эксплуатации скважин, геофизических исследований скважин» [32]. Уравнения, описывающие динамику пластовой системы, в общем случае являются многомерными нестационарными нелинейными уравнениями в частных производных. Возникающие краевые задачи, как правило, очень сложны, решать их аналитическими методами затруднительно, поэтому обычно прибегают к математическому моделированию [65, 66]. Построение полноразмерной модели натыкается на недостаток достоверной информации о пластовой системе и происходящих в ней процессах для каждой конкретной добывающей скважины. 21

Поэтому от модели требуется отображения качественных закономерностей и таких количественных показателей, которые устойчивы к вариациям исходных данных. При это следует исходить из идеи «удовлетворительной» («разумной») точности моделирования. минимизировать Возможную включением в потерю модель точности адаптационных можно механизмов, обеспечивающих гибкую настройку и изменение при необходимости модели под каждый конкретный рабочий режим скважины и изменившийся характер нефтедобычи. К математической модели «скважина - нефтяной пласт - скважина» помимо отражения реальных физических процессов, происходящих в нефтяном пласте и скважине при добыче, добавляются требования, связанные с технологией математического моделирования. Математическая модель должна просто настраиваться и изменяться при возникновении изменений в условиях добычи, параметрах пласта и скважины. Для математического моделирования работы системы «скважина – пласт – скважина» [57, 58, 77, 78] успешно может применяться технология емкостных моделей (Capacitance Model) [34, 70-75], которая представляет собой модель «вход-выход», характеризующая свойства нефтегазвого пласта исключительно на основании использования данных истории по динамике добычи и нагнетания. В емкостной модели сигнал входа (интенсивность нагнетания) преобразуется в сигнал выхода (дебит). Форма ответа выхода на добывающей скважине при ступенчатом изменении скорости нагнетания на входе в нагнетательной скважине зависит от временных задержек и затухания сигнала, которые зависят от среды между добывающей и нагнетательной скважинами. Предложенная технология [77] позволяет измерить межскважинную проводимость и степень накопления жидкости. Поскольку для нахождения параметров емкостной модели и оценки связности скважин требуются только данные по нагнетанию жидкости, дебиту скважин, которые обычно уже измерены и известны, то эта модель предлагают проводить быструю оценку связности между нагнетательными и добывающими скважинами по сравнению с традиционными моделями пластов. Также описываемый метод не требует 22

предварительного построения геологической модели исследуемого нефтяного месторождения. На сегодняшний день существует множество подходов и методик для такого рода анализа. Здесь можно выделить подход с расчетом коэффициентов корреляции между скоростями добычи и закачки на соседних скважинах [36]. Предложенный подход довольно прост в использовании, однако не учитывает интерференцию влияний нескольких соседних нагнетательных скважин на одну и ту же добывающую скважину. Также существует метод анализа, основанный на матрице взаимной продуктивности, который был предложен [35]. Его особенностью является сильная зависимость от качества исходных данных по забойным давлениям на скважинах. Наиболее точные и надежные результаты в этом методе получаются в случае, если скважины часто отключались в процессе эксплуатации на временном интервале по которому используются данные эксплуатации. В случае если забойные давление на всем используемом временном интервале не изменялись и скважины не отключались, то решение становится неустойчивым. В работе [33] предложен улучшенный метод расчета межскважинной связности пласта. Описываемый метод основан на объединении двух подходов – матрицы взаимной продуктивности (MPI – Multiwell Productivity Index) и емкостной модели (CM – Capacitance Model). Здесь MPI рассчитывается на основе геометрии пласта и скважин с усредненными свойствами, без использования данных эксплуатации. Анализ связности скважин производится в рамках емкостной модели (Capacitance Model) многоскважинной системы. Данная модель представляет дебиты на добывающих скважинах как функцию от состояния и истории их окружения, что позволяет количественно оценивать влияние соседних скважин друг на друга и предсказывать добычу по скважинам и месторождению в целом. Базовым уравнением для вывода уравнения емкостной модели многоскважинной системы является уравнение материального баланса для 23

системы пары скважин нагнетательная-добывающая (1.4.6) и уравнение Дюпюи (1.4.7) С𝑡𝑡 𝑉𝑉𝑝𝑝 𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑 = 𝑤𝑤(𝑡𝑡) − 𝑞𝑞(𝑡𝑡), 𝑞𝑞 = 𝐽𝐽(𝑝𝑝𝑞𝑞 − 𝑝𝑝̅), где: − С𝑡𝑡 – полная сжимаемость; − 𝑉𝑉𝑝𝑝 – дренируемый поровый объем; − 𝑝𝑝̅ – среднее пластовое давление; − 𝑝𝑝𝑞𝑞 – забойное давление добывающей скважины; − 𝐽𝐽 – индекс продуктивности добывающей скважины; − 𝑞𝑞(𝑡𝑡), 𝑤𝑤(𝑡𝑡) – скорости добычи и нагнетания жидкости. В работе [37] основная развиваемая идея заключается в использовании раннее упомянутой матрицы взаимной продуктивности [35]. Что позволяет связать вектор дебита скважин с вектором депрессий. Для случая, когда пласт вскрыт только одной скважиной, распределение давления в пласте при псевдоустановившемся режиме было получено [38]. Решение для системы из нескольких добывающих и нагнетательных скважин получается с применением принципа суперпозиции, принимая во внимание перепад давления, обусловленный наличием скин-фактора. Элементы вектора скин-фактора, отражают не только механический скин-фактор данной скважины, а скорее являются мерой отклонения производительности скважины от идеальной. Имея данные по истории добычи для группы скважин (временные ряды давлений и дебитов по всем рассматриваемым скважинам) можно реализовать сравнительный анализ производительности каждой скважины относительно других скважин. Описываемая идея применима только в многоскважинных системах. В этом случае можно представить, что явления, связанные с пластом, влияют на все скважины в равной мере. Сравнение графиков динамики скинфактора позволяет выявить проблемы отдельных скважин, поведение которых 24

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

ГЛАВА 2. ПАССИВНАЯ ГИДРОДИНАМИЧЕСКАЯ ТОМОГРАФИЯ Глава содержит описание построения математической модели эксплуатации нефтегазового месторождения на основе технологии емкостного моделирования с использованием данных штатной эксплуатации скважин месторождения для синтеза данных временных интервалов восстановления давления между реконструкция парами скважинами, коэффициента по которым пьезопроводности осуществляется проницаемого пласта. Предложенная модель [103, 104] сможет синтезировать необходимые исходные данные для гидродинамической томографии, путем имитации депрессии в скважинах и веерную регистрацию откликов на созданную депрессию с соседних скважин. 2.1. Математическая модель динамики эксплуатации месторождения Построение модели динамики работы нефтегазового месторождения для ее последующего использования с целью синтеза интервальных времен 𝜏𝜏𝑖𝑖,𝑗𝑗 на основе выполнения томографии (вместо вычислительного натурного) эксперимента основано на гидродинамической технологии емкостной модели (Capacitance Model). Выбранная технология моделирования характеризует свойства продуктивного пласта с использованием данных из истории штатной эксплуатации нефтегазового месторождения: объемы (или скорости) по дебиту (скорость добычи углеводородов на скважине в единицу времени) и нагнетания (закачка в скважину) жидкости. Построена модель гидродинамической связи нефтегазового месторождения, в котором происходит закачка и добыча углеводородов на интервалах времени, адаптированной к поставленной задаче гидродинамической томографии. Математическая модель месторождения представляет собой системы связанных скважин, где дебит скважин интервале времени находятся с применением принципа суперпозиции трех отдельных физических процессов, 26

проходящих внутри проницаемого пласта. При построении модели было принято решение руководствоваться простыми базовыми принципами и хорошими аппроксимационными способностями построенной конструкции модели [47]: 𝑄𝑄𝑖𝑖 (𝑡𝑡) = 𝑄𝑄1 (𝑡𝑡) + 𝑄𝑄2,𝑖𝑖 (𝑡𝑡) + 𝑄𝑄3,𝑖𝑖 (𝑡𝑡) Первый физический фактор 𝑄𝑄𝑖𝑖,1 (𝑡𝑡) это процесс, реализующий динамику начального дебита добывающей скважины i, который не подвержен влиянию работы соседних скважин, и аппроксимируется с применением линейного эволюционного уравнения: 𝑄𝑄1,𝑖𝑖 (𝑡𝑡) = 𝑒𝑒 −𝑡𝑡𝜆𝜆𝑖𝑖 𝑄𝑄𝑖𝑖,1 (𝑡𝑡0 ), здесь 𝑡𝑡0 – первоначальное время, когда запускается скважина на рассматриваемом интервале, 𝜆𝜆𝑖𝑖 – коэффициент затухания, от которого зависит скорость экспоненциального падения объем добычи углеводородов на скважине с номеров i и в частном случае зависит от времени, на котором определяется дебит. Следующий физический фактор 𝑄𝑄2,𝑖𝑖 (𝑡𝑡) это влияние окружающих нагнетательных скважин на дебит рассчитываемой i-ю добывающей скважины и рассчитывается как сумма влияний всех скважин, по которым происходит нагнетание жидкости в пласт: 𝑁𝑁 𝑖𝑖𝑖𝑖𝑖𝑖 𝑄𝑄𝑖𝑖,2 (𝑡𝑡) = ∑𝑗𝑗=1 Ψ𝑗𝑗 �𝑡𝑡 − 𝜎𝜎𝑖𝑖𝑖𝑖 �. Эта компонента модели определяет, как происходит перераспределение внутрипластового давления под влиянием проводимого нагнетания на месторождении, из-за чего изменяется фильтрационные течения внутри пласта и динамика перемещения углеводородов в связанной системе скважин. При определении влияния нагнетательных скважин стоит отметить, что только с определенной долей упрощения можно считать, что это влияние реализуется через линейную комбинацию притоков от нагнетательных скважин (2.1.5). Причем на величину притока от отдельной нагнетательной скважины 27

накладывает свое влияние коэффициент, учитывающий экспоненциальную задержку воздействия работы нагнетательной скважины во времени. Влияние отдельной j-й нагнетательной скважины находится по следующей формуле: Ψ𝑗𝑗 �𝑡𝑡, 𝜎𝜎𝑖𝑖𝑖𝑖 � = 𝛽𝛽𝑖𝑖𝑖𝑖 Ι𝑗𝑗 �𝑡𝑡 − 𝜎𝜎𝑖𝑖𝑖𝑖 �Ε𝑖𝑖𝑖𝑖𝑖𝑖 , 𝜎𝜎𝑖𝑖𝑖𝑖 = 𝑅𝑅𝑖𝑖𝑖𝑖 . 𝑉𝑉𝑖𝑖𝑖𝑖 Коэффициент запаздывания сигнала 𝜎𝜎𝑖𝑖𝑖𝑖 (2.1.7) зависит от расстояния на сетке между парой (𝑖𝑖, 𝑗𝑗) скважин – 𝑅𝑅𝑖𝑖𝑖𝑖 и скорости движения флюидов. Коэффициент 𝛽𝛽𝑖𝑖𝑗𝑗 часть интерференции (воздействия) j-й нагнетательной скважины на добывающую i-ю, с учетом того что суммарно коэффициенты интерференции по скважине j не могут превышать единицы по всем добывающим скважинам. Элемент 𝛪𝛪𝑗𝑗 – это дополнительное поступление флюидов к i-й добывающей скважине за счет работы j-я нагнетательной скважины к времени t: 𝑡𝑡 Ι𝑗𝑗 �𝑡𝑡 − 𝜎𝜎𝑖𝑖𝑖𝑖 � = � [𝑊𝑊𝑗𝑗 �𝑡𝑡𝑤𝑤 − 𝜎𝜎𝑖𝑖𝑖𝑖 � ∆𝑡𝑡] − Η𝑗𝑗 �𝑡𝑡 − 𝜎𝜎𝑖𝑖𝑖𝑖 �, 𝑡𝑡𝑤𝑤 =𝑡𝑡0 𝑡𝑡−1 𝑁𝑁𝑜𝑜𝑜𝑜𝑜𝑜 Η𝑗𝑗 (𝑡𝑡) = � � Ψ𝑗𝑗 �𝑡𝑡𝐻𝐻 − 𝜎𝜎𝑖𝑖𝑖𝑖 �. 𝑡𝑡𝐻𝐻 =𝑡𝑡0 𝑖𝑖=1 В математической модели 𝑁𝑁𝑖𝑖𝑖𝑖𝑖𝑖 и 𝑁𝑁𝑜𝑜𝑜𝑜𝑜𝑜 – количество добывающих и нагнетательных скважин, 𝛼𝛼𝑖𝑖𝑖𝑖𝑖𝑖 – коэффициент снижения (затухания) скорости движения флюидов, который представляет собой гидродинамическое сопротивление проницаемой среды перемещению притока от закачиваемой жидкости к участку добывающей скважины. На значение создаваемого притока (2.1.6) накладывает свое влияние скорость нагнетания жидкости 𝑊𝑊𝑗𝑗 в скважины и коэффициенты задержки сигнала (временной лаг) 𝜎𝜎𝑖𝑖𝑖𝑖 (2.1.7) между парой скважин. При этом в определении элемента Ι𝑗𝑗 не рассматривается влияние j-й скважины 𝐻𝐻𝑗𝑗 (2.1.7), которое уже был учтено при расчете дебита по 28

всем рассматриваемым добывающим скважинам к определенному моменту времени 𝑡𝑡. Элемент модели Ε𝑖𝑖𝑖𝑖𝑖𝑖 – это экспоненциальное затухание перемещения флюидов в проницаемом пласте, связанное с коэффициентом временного лага и реализуется по экспоненциальному закону: Ε𝑖𝑖𝑖𝑖𝑖𝑖 = 𝑒𝑒 −𝛼𝛼𝑖𝑖𝑖𝑖𝑖𝑖𝜎𝜎𝑖𝑖𝑖𝑖 Третий фактор дебита 𝑄𝑄3,𝑖𝑖 (𝑡𝑡) отвечает за влияние добычи углеводородов в соседних скважинах и также определяется линейной комбинацией по окружающим скважинам: 𝑄𝑄3,𝑖𝑖 (𝑡𝑡) = 𝑁𝑁𝑜𝑜𝑜𝑜𝑜𝑜 � 𝑗𝑗=1(𝑗𝑗≠𝑖𝑖) Φ𝑗𝑗 �𝑡𝑡 − 𝜎𝜎𝑖𝑖𝑖𝑖 � Значение Φ𝑗𝑗 выражает влияние дебита j-й добывающей скважины на дебит скважины i, рассчитывающееся по формуле: Φ𝑗𝑗 �𝑡𝑡, 𝜎𝜎𝑖𝑖𝑖𝑖 � = 𝛾𝛾𝑖𝑖𝑖𝑖 Δ𝐺𝐺𝑖𝑖𝑖𝑖 �𝑡𝑡 − 𝜎𝜎𝑖𝑖𝑖𝑖 �Ε𝑜𝑜𝑜𝑜𝑜𝑜 , Δ𝐺𝐺𝑖𝑖𝑖𝑖 = 𝐺𝐺𝑖𝑖 �𝑡𝑡 − 𝜎𝜎𝑖𝑖𝑖𝑖 � − 𝐺𝐺𝑗𝑗 �𝑡𝑡 − 𝜎𝜎𝑖𝑖𝑖𝑖 �, Ε𝑜𝑜𝑜𝑜𝑜𝑜 = 𝑒𝑒 −𝛼𝛼𝑜𝑜𝑜𝑜𝑜𝑜𝜎𝜎𝑖𝑖𝑖𝑖 . Где Δ𝐺𝐺𝑖𝑖𝑖𝑖 – разница в объёмах (или скорости) извлекаемых флюидов в скважинах i и j (2.1.12), 𝛾𝛾𝑖𝑖𝑖𝑖 – коэффициент влияния дебита скважины j на работу i-й скважины. Снижение скорости движения потоков флюидов Ε𝑜𝑜𝑜𝑜𝑜𝑜 (2.1.13) вычисляется через коэффициент задержки сигнала между двумя скважинами 𝜎𝜎𝑖𝑖𝑖𝑖 и параметр гидродинамического сопротивления перемещению обратного потока флюидов 𝛼𝛼𝑜𝑜𝑜𝑜𝑜𝑜 э. Общая аналитическая модель принимает следующую форму: 𝑁𝑁 𝑖𝑖𝑖𝑖𝑖𝑖 𝑄𝑄𝑖𝑖 (𝑡𝑡) = 𝑒𝑒 −𝑡𝑡𝜆𝜆𝑖𝑖 𝑄𝑄𝑖𝑖,1 (𝑡𝑡0 ) + ∑𝑗𝑗=1 [𝛽𝛽𝑖𝑖𝑖𝑖 [∑𝑡𝑡𝑡𝑡𝑤𝑤=𝑡𝑡0�𝑊𝑊𝑗𝑗 �𝑡𝑡𝑤𝑤 − 𝜎𝜎𝑖𝑖𝑖𝑖 �∆𝑡𝑡� − 𝑁𝑁 𝑜𝑜𝑜𝑜𝑜𝑜 Η𝑗𝑗 �t − 𝜎𝜎𝑖𝑖𝑖𝑖 �]𝑒𝑒 −𝛼𝛼𝑖𝑖𝑖𝑖𝑖𝑖𝜎𝜎𝑖𝑖𝑖𝑖 ] + ∑𝑗𝑗=1(𝑗𝑗≠𝑖𝑖) [𝛾𝛾𝑖𝑖𝑖𝑖 �𝐺𝐺𝑖𝑖 �𝑡𝑡 − 𝜎𝜎𝑖𝑖𝑖𝑖 � − 𝐺𝐺𝑗𝑗 �𝑡𝑡 − 𝜎𝜎𝑖𝑖𝑖𝑖 ��𝑒𝑒 −𝛼𝛼𝑜𝑜𝑜𝑜𝑜𝑜 𝜎𝜎𝑖𝑖𝑖𝑖 ]. 29

Представим уравнение (14) в символьной форме: 𝑄𝑄𝑖𝑖 (𝑡𝑡) = 𝐴𝐴[𝛼𝛼𝑖𝑖𝑖𝑖𝑖𝑖 , 𝛼𝛼𝑜𝑜𝑜𝑜𝑜𝑜 , 𝜆𝜆, 𝛽𝛽, 𝛾𝛾, 𝑉𝑉] Построение математической модели динамики работы нефтегазового месторождения реализуется через прогноз работы скважин по истории их эксплуатации, используя следующие параметры для оптимизации модели: 𝛼𝛼𝑖𝑖𝑖𝑖𝑖𝑖 , 𝛼𝛼𝑜𝑜𝑜𝑜𝑜𝑜 , 𝜆𝜆 = {𝜆𝜆𝑖𝑖 }, 𝛾𝛾 = �𝛾𝛾𝑖𝑖𝑖𝑖 �, 𝛽𝛽 = �𝛽𝛽𝑖𝑖𝑖𝑖 �, 𝑉𝑉 = �𝑉𝑉𝑖𝑖𝑖𝑖 �. Исходными данными для оптимизации (настройки модели) являются значения: − {𝐺𝐺𝑖𝑖 (𝑡𝑡)}; − история разработки: нагнетание - 𝑊𝑊 = �𝑊𝑊𝑗𝑗 (𝑡𝑡)�, добыча - 𝐺𝐺 = нагнетательные и добывающие скважины 𝑁𝑁𝑜𝑜𝑜𝑜𝑜𝑜 и 𝑁𝑁𝑖𝑖𝑖𝑖𝑖𝑖 . Модельные параметров среды выбираются для каждой рассматриваемой системы скважин нефтегазового месторождения так, чтобы уже известная история штатной эксплуатации скважин при подстановке в модель (2.1.14) давала историю дебитов с минимальной погрешностью относительно реальной истории эксплуатации. Следовательно, для расчета подходящих модельных параметров необходимо решить следующую оптимизационную задачу: 𝑇𝑇 𝑁𝑁𝑜𝑜𝑜𝑜𝑜𝑜 𝑍𝑍(𝛼𝛼𝑖𝑖𝑖𝑖𝑖𝑖 , 𝛼𝛼𝑜𝑜𝑜𝑜𝑜𝑜 , 𝜆𝜆, 𝛽𝛽, 𝛾𝛾, 𝑉𝑉) = � � |𝑄𝑄�𝚤𝚤 (𝑡𝑡) − 𝑄𝑄𝑖𝑖 (𝑡𝑡)| → 𝑚𝑚𝑚𝑚𝑚𝑚 𝑡𝑡=𝑡𝑡0 𝑖𝑖=1 Здесь 𝑄𝑄�𝑖𝑖 – дебит скважины 𝑖𝑖 из истории эксплуатации скважин и 𝑄𝑄𝑖𝑖 - модельный дебит скважины. Решение модели реализуется с учетом определенных условий накладываемый на некоторые параметры модели, связанных с физическим смыслом описываемых явлений. Условие по ограничению параметра 𝛽𝛽𝑖𝑖𝑗𝑗 зависит от того, что от скважины, нагнетающей жидкость в пласт, должно поступать притока в добывающие скважины не больше, чем было закачано. Поскольку 𝛽𝛽𝑖𝑖𝑖𝑖 определяет, какая часть жидкости из j-й нагнетательной 30

скважины движется в сторону добывающей скважины i, то он должен удовлетворять следующему ограничению: 𝑁𝑁𝑜𝑜𝑜𝑜𝑜𝑜 � 𝛽𝛽𝑖𝑖𝑗𝑗 ≤ 1, для всех 𝑗𝑗 = 1 … 𝑁𝑁𝑖𝑖𝑖𝑖𝑖𝑖 𝑖𝑖=1 Величина параметра определяется в границе от 0 до 1, как доля жидкости, нагнетаемой в скважину 𝑗𝑗 и влияющей на работу добывающей скважины. 0 ≤ 𝛽𝛽𝑖𝑖𝑖𝑖 ≤ 1 Ограничение параметра 𝛾𝛾𝑖𝑖𝑖𝑖 обусловлено смыслом третьей компоненты модели, как конкуренции добывающих скважин за приток флюидов внутри пласта и не может привносить дополнительную массу в систему скважин. Поэтому сумма 𝛾𝛾𝑖𝑖𝑖𝑖 по всем скважинам должна быть равно нулю в рассматриваемый момент времени. 𝑁𝑁𝑜𝑜𝑜𝑜𝑜𝑜 𝑁𝑁𝑜𝑜𝑜𝑜𝑜𝑜 � � 𝛾𝛾𝑖𝑖𝑗𝑗 = 0 𝑖𝑖(𝑗𝑗)=1 𝑗𝑗(𝑖𝑖)=1 Параметр определяет степень влияния между добывающими скважинами, то значения его должны лежать в интервале от 0 до 1: Параметры калибровочными 0 ≤ 𝛾𝛾𝑖𝑖𝑖𝑖 ≤ 1 математической коэффициентами, модели от 𝛼𝛼𝑖𝑖𝑖𝑖𝑖𝑖 , 𝛼𝛼𝑜𝑜𝑜𝑜𝑜𝑜 , 𝜆𝜆𝑖𝑖𝑖𝑖 которых зависит являются скорость экспоненциального замедления протекания физических процессов внутри проницаемого пласта, которые они определяют, а их величины находятся в границах от 0 (нет влияния на процесс) до 1. Проверка на искусственном примере. Разработанная математическая модель многоскважинной системы была проверена на искусственном примере разработки месторождения. Была выбрана сетка скважин с несколькими скважинами: четыре нагнетательные скважины расположены по углам сетки и четыре добывающих между ними (Рисунок 2.1.1). Продолжительность разработки месторождения составила 30 месяцев (Рисунок 2.1.2). По тестовым 31

данным истории разработки построена модель месторождения. В таком численном эксперименте проверялась способность модели подбить модельные дебиты к дебитам из тестовой эксплуатации скважин нефтегазового месторождения (Рисунки 2.1.3, 2.1.4, 2.1.5). Сетка скважин тестового месторождения История разработки месторождения (добыча – верхняя полуплоскость, нагнетание – нижняя полуплоскость) 32

а) Нулевое приближение перед оптимизацией б) После оптимизации Динамика реального и модельного дебитов на рассматриваемом интервале времени 33

а) Нулевое приближение перед б) После оптимизации оптимизацией Динамика реального и модельного дебитов на рассматриваемом интервале времени для скважины №1 а) Нулевое приближение перед б) После оптимизации оптимизацией Динамика реального и модельного дебитов на рассматриваемом интервале времени для скважины №2 Модельные дебиты достаточно близко переместились к значениям дебитов по скважинам из истории экспериментального месторождения, и относительная ошибка по добывающим скважинами составила 5, 87 %. Полученный результат оптимизации 34 модели можно считать

удовлетворительным, то есть на основе полученных параметров взаимодействия между скважинами можно будет сделать прогноз динамики добычи жидкости по каждой скважине. Так что построенную математическую модель можно использовать для синтезирования начальных данных в виде измерения временных интервалов наступления реакции в скважинам- приемниках на изменение режима работы в скважинах-источниках для метода гидродинамической томографии. 2.2. Итерационный алгоритм подбора оптимального вектора параметров По имеющимся данным из истории эксплуатации месторождения, представляющей собой временной поток данных (дебитов: положительных – работа добывающих скважин, отрицательных – работа нагнетательных скважин) системы из N скважин построить ее математическую модель. 1) 𝑦𝑦� – наблюдаемые данные. 𝐴𝐴�𝑋𝑋�, 𝜉𝜉 ̅� = 𝑦𝑦� (2.2.1) 2) А – оператор прогноза по математической модели эксплуатации месторождения. 3) 𝑋𝑋� – вектор входных данных - – история штатной эксплуатации месторождения. 4) 𝜉𝜉 ̅ – вектор параметров оптимизации – вектор параметров по которому осуществляется оптимизация. В качестве входных данных для модели прогноза служат 𝑋𝑋� и 𝜉𝜉 ̅: 1. На временном ряде от начального интервала 𝑡𝑡0 до конечного интервала 𝑡𝑡𝑛𝑛 заданы значения дебитов работы скважин в виде временного ряда: 𝑋𝑋� = �𝐺𝐺𝑖𝑖 (𝑡𝑡), 𝑊𝑊𝑗𝑗 (𝑡𝑡)�. 𝑖𝑖 = 1 … 𝑁𝑁𝑖𝑖𝑖𝑖𝑖𝑖 , 𝑗𝑗 = 1 … 𝑁𝑁𝑜𝑜𝑜𝑜𝑜𝑜 , 𝑡𝑡 = 𝑡𝑡0 … 𝑡𝑡𝑛𝑛 a. 𝑛𝑛 – количество временных интервалов истории разработки. b. 𝐺𝐺𝑖𝑖 (𝑡𝑡) – дебит добывающей скважины i на интервале времени 𝑡𝑡ℎ . 𝑁𝑁𝑖𝑖𝑖𝑖𝑖𝑖 – количество добывающих скважин. 35

c. 𝑊𝑊𝑗𝑗 (𝑡𝑡) - закачка нагнетательной скважины j на интервале времени 𝑡𝑡ℎ . 𝑁𝑁𝑜𝑜𝑜𝑜𝑜𝑜 – количество нагнетательных скважин. 2. Вектор оптимизационных параметров - 𝜉𝜉 ̅ = {𝛼𝛼𝑖𝑖𝑖𝑖𝑖𝑖 , 𝛼𝛼𝑜𝑜𝑜𝑜𝑜𝑜 , 𝜆𝜆, 𝛽𝛽, 𝛾𝛾, 𝑉𝑉}. Решение задачи состоит в подборе вектора параметров модели 𝜉𝜉 ̅𝑛𝑛,𝑘𝑘 при построении оператора A(1) по выборкам из всех либо части предыдущих элементов из временного ряда данных истории разработки: �𝐺𝐺𝑖𝑖 (𝑡𝑡𝑚𝑚 ), 𝑊𝑊𝑗𝑗 (𝑡𝑡𝑚𝑚 )�, 𝑚𝑚 = 𝑛𝑛 − 𝑘𝑘 ÷ 𝑛𝑛, которые играют роль входного вектора параметров 𝑋𝑋�(𝑡𝑡𝑛𝑛 , 𝑘𝑘). k – количество временных интервалов используемых в прогнозе (k = 1 ÷ 𝑛𝑛) Процесс оптимизации модели состоит в подборе параметра 𝜉𝜉 ̅𝑛𝑛,𝑘𝑘 для минимизации невязки: ��𝑦𝑦(𝑡𝑡𝑛𝑛+1 ) − 𝐴𝐴[𝑋𝑋�(𝑡𝑡𝑛𝑛+1 , 𝑘𝑘), 𝜉𝜉 ̅𝑛𝑛,𝑘𝑘 ]�� → 𝑚𝑚𝑚𝑚𝑚𝑚 (2.2.2) 𝐴𝐴�𝑋𝑋�(𝑡𝑡𝑛𝑛+1 , 𝑘𝑘), 𝜉𝜉 ̅𝑛𝑛,𝑘𝑘 � = 𝑄𝑄�𝛿𝛿 (𝑡𝑡𝑛𝑛+1 , 𝑘𝑘) (2.2.3) Оператор прогноза: Уравнение (2.2.3) – результат реального вычисления по (2.2.1) с подобранными исходя из (2.2.2) оптимизационными параметрами 𝜉𝜉 ̅𝑛𝑛,𝑘𝑘 для выборки из истории разработки 𝑋𝑋� = {𝐺𝐺(𝑡𝑡𝑛𝑛 , 𝑘𝑘), 𝑊𝑊(𝑡𝑡𝑛𝑛+1 , 𝑘𝑘)}. соответствует 𝑦𝑦(𝑡𝑡𝑛𝑛+1 ) с уровнем погрешности 𝛿𝛿(𝑛𝑛, 𝑘𝑘): 𝛿𝛿(𝑛𝑛, 𝑘𝑘) = ||𝑄𝑄�𝛿𝛿 (𝑡𝑡𝑛𝑛+1 , 𝑘𝑘) − 𝑦𝑦�(𝑡𝑡𝑛𝑛+1 )||𝑅𝑅𝑁𝑁 𝑄𝑄�𝛿𝛿 (𝑡𝑡𝑛𝑛+1 , 𝑘𝑘) (2.2.4) Величина 𝛿𝛿(𝑛𝑛, 𝑘𝑘) характеризует предельно достигаемую динамику ошибки при наращивании временного интервала данных (по параметру 𝒏𝒏) и по глубине одновременно используемых временных данных при прогнозе очередного шага временной последовательности (k). Следует построить и адаптировать модель (подобрать параметры 𝜉𝜉 ̅𝑛𝑛,𝑘𝑘 ) к реальным данным истории эксплуатации месторождения. Итерационный метод построения решения для задачи [54, 55]: � 𝐴𝐴�𝑋𝑋�(𝑡𝑡𝑛𝑛+1 , 𝑘𝑘), 𝜉𝜉 ̅𝑛𝑛,𝑘𝑘 � = 𝑄𝑄�𝛿𝛿 (𝑡𝑡𝑛𝑛+1 , 𝑘𝑘) ��𝑦𝑦(𝑡𝑡𝑛𝑛+1 ) − 𝑄𝑄�𝛿𝛿 (𝑡𝑡𝑛𝑛+1 , 𝑘𝑘)�� 36 𝑅𝑅𝑁𝑁 → 𝑚𝑚𝑚𝑚𝑚𝑚 (2.2.5)

состоит в нахождении на очередной итерации вектора приращения по параметрам и проверки погрешности модели с новым вектором параметров. 𝜉𝜉 𝒏𝒏+𝟏𝟏 = 𝜉𝜉 𝑛𝑛 + 𝛼𝛼 𝑛𝑛 𝐴𝐴′∗ [𝑋𝑋, 𝜉𝜉 𝑛𝑛 ]𝜑𝜑 𝑛𝑛 𝜑𝜑𝑛𝑛 = 𝐴𝐴[𝑋𝑋, 𝜉𝜉 𝑛𝑛 ] − 𝑦𝑦. ‖𝐴𝐴′∗ [𝑋𝑋, 𝜉𝜉 𝑛𝑛 ]𝜑𝜑 𝑛𝑛 ‖2 𝛼𝛼 = ′ ‖𝐴𝐴 [𝑋𝑋, 𝜉𝜉 𝑛𝑛 ]𝐴𝐴′∗ [𝑋𝑋, 𝜉𝜉 𝑛𝑛 ]𝜑𝜑 𝑛𝑛 ‖2 𝑛𝑛 ′[ 𝐴𝐴 𝑋𝑋, 𝜉𝜉 𝑛𝑛 ] = Здесь: 𝐴𝐴′ [𝑋𝑋, 𝜉𝜉 𝑛𝑛 ]– матрица (2.2.6) 𝛿𝛿𝐴𝐴1 [𝑋𝑋,𝜉𝜉 𝑛𝑛 ] 𝛿𝛿𝐴𝐴1 [𝑋𝑋,𝜉𝜉 𝑛𝑛 ] 𝛿𝛿𝜉𝜉1𝑛𝑛 𝛿𝛿𝜉𝜉2𝑛𝑛 ⎧ 𝛿𝛿𝜉𝜉1𝑛𝑛 ⎪ 𝛿𝛿𝐴𝐴2[𝑋𝑋,𝜉𝜉𝑛𝑛 ] ⎨ … ⎪𝛿𝛿𝐴𝐴𝑀𝑀[𝑋𝑋,𝜉𝜉𝑛𝑛 ] ⎩ 𝛿𝛿𝜉𝜉1𝑛𝑛 {𝑴𝑴 × 𝑲𝑲} … 𝛿𝛿𝜉𝜉2𝑛𝑛 𝛿𝛿𝐴𝐴2 [𝑋𝑋,𝜉𝜉 𝑛𝑛 ] 𝛿𝛿𝐴𝐴 … 𝑀𝑀 [ 𝑋𝑋,𝜉𝜉 𝑛𝑛 ] 𝛿𝛿𝜉𝜉2𝑛𝑛 … … обратного 𝛿𝛿𝐴𝐴1 [𝑋𝑋,𝜉𝜉 𝑛𝑛 ] ⎫ ⎪ 𝛿𝛿𝐴𝐴2 [𝑋𝑋,𝜉𝜉 𝑛𝑛 ] 𝛿𝛿𝐴𝐴 𝛿𝛿𝜉𝜉𝐾𝐾𝑛𝑛 𝛿𝛿𝜉𝜉𝐾𝐾𝑛𝑛 … 𝑀𝑀 [ 𝑋𝑋,𝜉𝜉 𝑛𝑛 ] 𝛿𝛿𝜉𝜉𝐾𝐾𝑛𝑛 оператора производных по параметрам 𝜉𝜉 = {𝛼𝛼𝑖𝑖𝑖𝑖𝑖𝑖 , 𝛼𝛼𝑜𝑜𝑜𝑜𝑜𝑜 , 𝜆𝜆, 𝛽𝛽, 𝛾𝛾, 𝑉𝑉1 , 𝑉𝑉2 }. ⎬ ⎪ ⎭ А (2.2.7) в частных 𝐴𝐴′∗ [𝑋𝑋, 𝜉𝜉 𝑛𝑛 ]- сопряженная матрица {𝑲𝑲 × 𝑴𝑴} Если месторождение состоит из 4 скважин (2 – добывающие, 2 – нагнетательные), то матрица будет иметь размерность ({2 × 16}). 2 – количество добывающих скважин (M). 16 – количество параметров среды (K). Частные производные оператора A по параметрам 𝜉𝜉 = {𝛼𝛼𝑖𝑖𝑖𝑖𝑖𝑖 , 𝛼𝛼𝑜𝑜𝑜𝑜𝑜𝑜 , 𝜆𝜆, 𝛽𝛽, 𝛾𝛾, 𝑉𝑉1 , 𝑉𝑉2 } для итерационного процесса. Аналитическое решение в частных производных для оператора 𝐴𝐴, который имеет вид: 𝑁𝑁𝑖𝑖𝑖𝑖𝑖𝑖 𝑡𝑡 𝑗𝑗=1 𝑡𝑡𝑤𝑤 =𝑡𝑡0 𝐴𝐴[𝑋𝑋, 𝜉𝜉] = 𝑒𝑒 −𝑡𝑡𝜆𝜆𝑖𝑖 𝑄𝑄𝑖𝑖,1 (𝑡𝑡0 ) 𝑅𝑅𝑖𝑖,𝑗𝑗 𝑅𝑅𝑖𝑖,𝑗𝑗 𝑅𝑅𝑖𝑖,𝑗𝑗 −𝛼𝛼𝑖𝑖𝑖𝑖𝑖𝑖 𝑉𝑉𝑖𝑖,𝑗𝑗 ] + �[𝛽𝛽𝑖𝑖,𝑗𝑗 ( � [𝑊𝑊𝑗𝑗 �𝑡𝑡𝑤𝑤 − � ∆𝑡𝑡] − Η𝑗𝑗 �𝑡𝑡 − �)𝑒𝑒 𝑉𝑉𝑖𝑖,𝑗𝑗 𝑉𝑉𝑖𝑖,𝑗𝑗 + 𝑁𝑁𝑜𝑜𝑜𝑜𝑜𝑜 � 𝑖𝑖2=1(𝑖𝑖2≠𝑖𝑖) 𝑅𝑅 𝑅𝑅𝑖𝑖,𝑖𝑖2 𝑅𝑅𝑖𝑖,𝑖𝑖2 −𝛼𝛼𝑜𝑜𝑜𝑜𝑜𝑜 𝑖𝑖,𝑖𝑖2 𝑉𝑉𝑖𝑖,𝑖𝑖2 � � − 𝐺𝐺𝑖𝑖2 �𝑡𝑡 − �]𝑒𝑒 �𝛾𝛾𝑖𝑖,𝑖𝑖2 [𝐺𝐺𝑖𝑖 �𝑡𝑡 − 𝑉𝑉𝑖𝑖,𝑖𝑖2 𝑉𝑉𝑖𝑖,𝑖𝑖2 1. Производная по параметру 𝛼𝛼𝑖𝑖𝑖𝑖𝑖𝑖 : 37

𝛿𝛿𝛿𝛿[𝑋𝑋, 𝜉𝜉] = (𝑄𝑄1 + 𝑄𝑄2 + 𝑄𝑄3 )′ 𝛿𝛿𝛼𝛼𝑖𝑖𝑖𝑖𝑖𝑖 𝑁𝑁𝑖𝑖𝑖𝑖𝑖𝑖 𝑡𝑡 = 0 + ��[𝛽𝛽𝑖𝑖,𝑗𝑗 ( � [𝑊𝑊𝑗𝑗 �𝑡𝑡𝑤𝑤 − 𝑗𝑗=1 𝑁𝑁𝑖𝑖𝑖𝑖𝑖𝑖 𝑡𝑡𝑤𝑤=𝑡𝑡0 𝑡𝑡 𝑅𝑅𝑖𝑖,𝑗𝑗 𝑅𝑅𝑖𝑖,𝑗𝑗 � ∆𝑡𝑡] − Η𝑗𝑗 �𝑡𝑡 − �)𝑒𝑒 𝑉𝑉𝑖𝑖,𝑗𝑗 𝑉𝑉𝑖𝑖,𝑗𝑗 𝑅𝑅𝑖𝑖,𝑗𝑗 −𝛼𝛼𝑖𝑖𝑖𝑖𝑖𝑖 𝑉𝑉𝑖𝑖,𝑗𝑗 ′ ]� + 0 𝑅𝑅𝑖𝑖,𝑗𝑗 −𝛼𝛼𝑖𝑖𝑖𝑖𝑖𝑖𝑅𝑅𝑉𝑉𝑖𝑖,𝑗𝑗 𝑅𝑅𝑖𝑖,𝑗𝑗 𝑅𝑅𝑖𝑖,𝑗𝑗 𝑖𝑖,𝑗𝑗 ] )𝑒𝑒 = �[𝛽𝛽𝑖𝑖,𝑗𝑗 ( � [𝑊𝑊𝑗𝑗 �𝑡𝑡𝑤𝑤 − � ∆𝑡𝑡] − Η𝑗𝑗 �𝑡𝑡 − �)(− 𝑉𝑉𝑖𝑖,𝑗𝑗 𝑉𝑉𝑖𝑖,𝑗𝑗 𝑉𝑉𝑖𝑖,𝑗𝑗 𝑗𝑗=1 𝑡𝑡𝑤𝑤=𝑡𝑡0 2. Производная по параметру 𝛼𝛼𝑜𝑜𝑜𝑜𝑜𝑜 : 𝛿𝛿𝛿𝛿[𝑋𝑋, 𝜉𝜉] = (𝑄𝑄1 + 𝑄𝑄2 + 𝑄𝑄3 )′ 𝛿𝛿𝛼𝛼𝑜𝑜𝑜𝑜𝑜𝑜 =0+0+� = 𝑁𝑁𝑜𝑜𝑜𝑜𝑜𝑜 � 𝑖𝑖2=1(𝑖𝑖2≠𝑖𝑖) 𝑁𝑁𝑜𝑜𝑜𝑜𝑜𝑜 � 𝑖𝑖2=1(𝑖𝑖2≠𝑖𝑖) �𝛾𝛾𝑖𝑖,𝑖𝑖2 [𝐺𝐺𝑖𝑖 �𝑡𝑡 − 𝑅𝑅 𝑅𝑅𝑖𝑖,𝑖𝑖2 −𝛼𝛼𝑜𝑜𝑜𝑜𝑜𝑜 𝑖𝑖,𝑖𝑖2 𝑉𝑉𝑖𝑖,𝑖𝑖2 �� ]𝑒𝑒 𝑅𝑅𝑖𝑖,𝑖𝑖2 � − 𝐺𝐺𝑖𝑖2 �𝑡𝑡 − 𝑉𝑉𝑖𝑖,𝑖𝑖2 𝑉𝑉𝑖𝑖,𝑖𝑖2 𝑅𝑅𝑖𝑖,𝑖𝑖2 𝑅𝑅𝑖𝑖,𝑖𝑖2 𝑅𝑅𝑖𝑖,𝑖𝑖2 −𝛼𝛼𝑜𝑜𝑜𝑜𝑜𝑜𝑅𝑅𝑉𝑉𝑖𝑖,𝑖𝑖2 𝑖𝑖,𝑖𝑖2 � )𝑒𝑒 � − 𝐺𝐺𝑖𝑖2 �𝑡𝑡 − �](− �𝛾𝛾𝑖𝑖,𝑖𝑖2 [𝐺𝐺𝑖𝑖 �𝑡𝑡 − 𝑉𝑉𝑖𝑖,𝑖𝑖2 𝑉𝑉𝑖𝑖,𝑖𝑖2 𝑉𝑉𝑖𝑖,𝑖𝑖2 3. Производная по параметру 𝜆𝜆: ′ 𝛿𝛿𝐴𝐴[𝑋𝑋, 𝜉𝜉] = (𝑄𝑄1 + 𝑄𝑄2 + 𝑄𝑄3 )′ = �𝑒𝑒 −𝑡𝑡𝜆𝜆𝑖𝑖 𝑄𝑄𝑖𝑖,1 (𝑡𝑡0 )� + 0 + 0 = −𝑡𝑡𝑒𝑒 −𝑡𝑡𝜆𝜆𝑖𝑖 𝑄𝑄𝑖𝑖,1 (𝑡𝑡0 ) 𝛿𝛿𝜆𝜆𝑖𝑖 4. Производная по параметру 𝛽𝛽: 𝛿𝛿𝛿𝛿[𝑋𝑋, 𝜉𝜉] = (𝑄𝑄1 + 𝑄𝑄2 + 𝑄𝑄3 )′ 𝛿𝛿𝛽𝛽𝑖𝑖𝑖𝑖 𝑁𝑁𝑖𝑖𝑖𝑖𝑖𝑖 𝑡𝑡 𝑗𝑗=1 𝑡𝑡𝑤𝑤=𝑡𝑡0 = 0 + ��[𝛽𝛽𝑖𝑖,𝑗𝑗 ( � [𝑊𝑊𝑗𝑗 �𝑡𝑡𝑤𝑤 − 𝑡𝑡 ′ 𝑅𝑅𝑖𝑖,𝑗𝑗 𝑅𝑅𝑖𝑖,𝑗𝑗 � ∆𝑡𝑡] − Η𝑗𝑗 �𝑡𝑡 − �)𝑒𝑒 𝑉𝑉𝑖𝑖,𝑗𝑗 𝑉𝑉𝑖𝑖,𝑗𝑗 𝑅𝑅𝑖𝑖,𝑗𝑗 𝑅𝑅𝑖𝑖,𝑗𝑗 𝑅𝑅𝑖𝑖,𝑗𝑗 −𝛼𝛼𝑖𝑖𝑖𝑖𝑖𝑖 𝑉𝑉𝑖𝑖,𝑗𝑗 = ( � [𝑊𝑊𝑗𝑗 �𝑡𝑡𝑤𝑤 − � ∆𝑡𝑡] − Η𝑗𝑗 �𝑡𝑡 − �)𝑒𝑒 𝑉𝑉𝑖𝑖,𝑗𝑗 𝑉𝑉𝑖𝑖,𝑗𝑗 𝑡𝑡𝑤𝑤=𝑡𝑡0 5. Производная по параметру 𝛾𝛾: 38 𝑅𝑅𝑖𝑖,𝑗𝑗 −𝛼𝛼𝑖𝑖𝑖𝑖𝑖𝑖 𝑉𝑉𝑖𝑖,𝑗𝑗 ′ ]� + 0

𝛿𝛿𝛿𝛿[𝑋𝑋, 𝜉𝜉] = (𝑄𝑄1 + 𝑄𝑄2 + 𝑄𝑄3 )′ 𝛿𝛿𝛾𝛾𝑖𝑖𝑖𝑖 =0+0 +� 𝑁𝑁𝑜𝑜𝑜𝑜𝑜𝑜 � 𝑖𝑖2=1(𝑖𝑖2≠𝑖𝑖) �𝛾𝛾𝑖𝑖,𝑖𝑖2 (𝐺𝐺𝑖𝑖 �𝑡𝑡 − 𝑅𝑅 𝑅𝑅𝑖𝑖,𝑖𝑖2 −𝛼𝛼𝑜𝑜𝑜𝑜𝑜𝑜 𝑖𝑖,𝑖𝑖2 𝑉𝑉𝑖𝑖,𝑖𝑖2 �� )𝑒𝑒 𝑅𝑅𝑖𝑖,𝑖𝑖2 � − 𝐺𝐺𝑖𝑖2 �𝑡𝑡 − 𝑉𝑉𝑖𝑖,𝑖𝑖2 𝑉𝑉𝑖𝑖,𝑖𝑖2 𝑅𝑅 𝑅𝑅𝑖𝑖,𝑖𝑖2 𝑅𝑅𝑖𝑖,𝑖𝑖2 −𝛼𝛼𝑜𝑜𝑜𝑜𝑜𝑜 𝑖𝑖,𝑖𝑖2 𝑉𝑉𝑖𝑖,𝑖𝑖2 = (𝐺𝐺𝑖𝑖 �𝑡𝑡 − � − 𝐺𝐺𝑖𝑖2 �𝑡𝑡 − �)𝑒𝑒 𝑉𝑉𝑖𝑖,𝑖𝑖2 𝑉𝑉𝑖𝑖,𝑖𝑖2 6. Производная по параметру 𝑉𝑉1 : 𝛿𝛿𝛿𝛿[𝑋𝑋, 𝜉𝜉] = (𝑄𝑄1 + 𝑄𝑄2 + 𝑄𝑄3 )′ 𝛿𝛿𝑉𝑉𝑖𝑖𝑖𝑖1 =0 𝑁𝑁𝑖𝑖𝑖𝑖𝑖𝑖 𝑡𝑡 𝑗𝑗=1 𝑡𝑡𝑤𝑤=𝑡𝑡0 + ��[𝛽𝛽𝑖𝑖,𝑗𝑗 ( � [𝑊𝑊𝑗𝑗 �𝑡𝑡𝑤𝑤 − +0 𝑡𝑡 𝑅𝑅𝑖𝑖,𝑗𝑗 𝑅𝑅𝑖𝑖,𝑗𝑗 � ∆𝑡𝑡] − Η𝑗𝑗 �𝑡𝑡 − �)𝑒𝑒 𝑉𝑉𝑖𝑖,𝑗𝑗 𝑉𝑉𝑖𝑖,𝑗𝑗 𝑅𝑅𝑖𝑖,𝑗𝑗 −𝛼𝛼𝑖𝑖𝑖𝑖𝑖𝑖 𝑉𝑉𝑖𝑖,𝑗𝑗 ′ ]� ′ 𝑅𝑅𝑖𝑖,𝑗𝑗 −𝛼𝛼𝑖𝑖𝑖𝑖𝑖𝑖𝑅𝑅𝑉𝑉𝑖𝑖,𝑗𝑗 𝑅𝑅𝑖𝑖,𝑗𝑗 𝑅𝑅𝑖𝑖,𝑗𝑗 𝑖𝑖,𝑗𝑗 = (𝛽𝛽𝑖𝑖,𝑗𝑗 ( � [𝑊𝑊𝑗𝑗 �𝑡𝑡𝑤𝑤 − )𝑒𝑒 � ∆𝑡𝑡] − Η𝑗𝑗 �𝑡𝑡 − �)(𝛼𝛼𝑖𝑖𝑖𝑖𝑖𝑖 𝑉𝑉𝑖𝑖,𝑗𝑗 𝑉𝑉𝑖𝑖,𝑗𝑗 𝑉𝑉𝑖𝑖,𝑗𝑗 2 𝑡𝑡𝑤𝑤=𝑡𝑡0 7. Производная по параметру 𝑉𝑉2 : =0+0 +� 𝑁𝑁𝑜𝑜𝑜𝑜𝑜𝑜 � 𝑖𝑖2=1(𝑖𝑖2≠𝑖𝑖) 𝛿𝛿𝛿𝛿[𝑋𝑋, 𝜉𝜉] = (𝑄𝑄1 + 𝑄𝑄2 + 𝑄𝑄3 )′ 𝛿𝛿𝑉𝑉𝑖𝑖𝑖𝑖2 �𝛾𝛾𝑖𝑖,𝑖𝑖2 (𝐺𝐺𝑖𝑖 �𝑡𝑡 − 𝑅𝑅 𝑅𝑅𝑖𝑖,𝑖𝑖2 −𝛼𝛼𝑜𝑜𝑜𝑜𝑜𝑜 𝑖𝑖,𝑖𝑖2 𝑉𝑉𝑖𝑖,𝑖𝑖2 �� )𝑒𝑒 𝑅𝑅𝑖𝑖,𝑖𝑖2 � − 𝐺𝐺𝑖𝑖2 �𝑡𝑡 − 𝑉𝑉𝑖𝑖,𝑖𝑖2 𝑉𝑉𝑖𝑖,𝑖𝑖2 ′ R i,i2 R i,i2 𝑅𝑅𝑖𝑖,𝑖𝑖2 −𝛼𝛼𝑜𝑜𝑜𝑜𝑜𝑜𝑅𝑅𝑉𝑉𝑖𝑖,𝑖𝑖2 𝑖𝑖,𝑖𝑖2 = 𝛾𝛾𝑖𝑖,𝑖𝑖2 (𝐺𝐺𝑖𝑖 �t − � − Gi2 �t − �)(𝛼𝛼𝑜𝑜𝑜𝑜𝑜𝑜 2 )𝑒𝑒 Vi,i2 Vi,i2 𝑉𝑉𝑖𝑖,𝑖𝑖2 Вектор параметров среды 𝜉𝜉 для модели из 4 скважин (2 добывающие скважины: №1 и №2, 2 нагнетательные скважины №3 и №4. 1 1 1 1 2 2 𝜉𝜉 = {𝛼𝛼𝑖𝑖𝑖𝑖𝑖𝑖 , 𝛼𝛼𝑜𝑜𝑜𝑜𝑜𝑜 , 𝜆𝜆1 , 𝜆𝜆2 , 𝛽𝛽13 , 𝛽𝛽14 , 𝛽𝛽23 , 𝛽𝛽24 , 𝛾𝛾12 , 𝛾𝛾21, 𝑉𝑉13 , 𝑉𝑉14 , 𝑉𝑉23 , 𝑉𝑉24 , 𝑉𝑉12 , 𝑉𝑉21 }. 𝐴𝐴′∗ [𝑋𝑋, 𝜉𝜉 𝑛𝑛 ] – матрица (16 × 2 (𝐾𝐾 × 𝑀𝑀)) – сопряженная к матрице обратного оператора 𝐴𝐴′ [𝑋𝑋, 𝜉𝜉 𝑛𝑛 ]: 39

𝛿𝛿𝐴𝐴1 [𝑋𝑋,𝜉𝜉 𝑛𝑛 ] 𝛿𝛿𝛼𝛼𝑖𝑖𝑖𝑖𝑖𝑖 ⎛𝛿𝛿𝐴𝐴1[𝑋𝑋,𝜉𝜉 𝑛𝑛] ⎜ 𝛿𝛿𝛼𝛼𝑜𝑜𝑜𝑜𝑜𝑜 ⎜𝛿𝛿𝐴𝐴1[𝑋𝑋,𝜉𝜉 𝑛𝑛] ⎜ 𝛿𝛿𝜆𝜆1 ⎜ 0 ⎜𝛿𝛿𝐴𝐴1[𝑋𝑋,𝜉𝜉 𝑛𝑛] ⎜ 𝛿𝛿𝛽𝛽13 ⎜𝛿𝛿𝐴𝐴1[𝑋𝑋,𝜉𝜉 𝑛𝑛] ⎜ 𝛿𝛿𝛽𝛽14 ⎜ 0 ′∗ [𝑋𝑋, 𝑛𝑛 ] 𝜉𝜉 = ⎜ 0 𝐴𝐴 ⎜𝛿𝛿𝐴𝐴1[𝑋𝑋,𝜉𝜉 𝑛𝑛] ⎜ 𝛿𝛿𝛾𝛾12 ⎜ 0 ⎜𝛿𝛿𝐴𝐴1[𝑋𝑋,𝜉𝜉 𝑛𝑛] ⎜ 𝛿𝛿𝑉𝑉13 ⎜𝛿𝛿𝐴𝐴1[𝑋𝑋,𝜉𝜉 𝑛𝑛] ⎜ 𝛿𝛿𝑉𝑉14 ⎜ 0 ⎜ 0 ⎜𝛿𝛿𝐴𝐴1[𝑋𝑋,𝜉𝜉 𝑛𝑛] ⎝ 𝛿𝛿𝑉𝑉12 0 𝛿𝛿𝐴𝐴2 [𝑋𝑋,𝜉𝜉 𝑛𝑛 ] 𝛿𝛿𝛼𝛼𝑖𝑖𝑖𝑖𝑖𝑖 𝛿𝛿𝐴𝐴2 [𝑋𝑋,𝜉𝜉 𝑛𝑛 ]⎞ ⎟ 0 ⎟ 2 𝛿𝛿𝐴𝐴 [𝑋𝑋,𝜉𝜉 𝑛𝑛 ]⎟ 𝛿𝛿𝜆𝜆2 ⎟ 0 ⎟ 0 ⎟ 2 𝛿𝛿𝐴𝐴 [𝑋𝑋,𝜉𝜉 𝑛𝑛 ]⎟ 𝛿𝛿𝛽𝛽23 ⎟ 𝛿𝛿𝐴𝐴2 [𝑋𝑋,𝜉𝜉 𝑛𝑛 ]⎟ 𝛿𝛿𝛽𝛽24 ⎟ 0 ⎟ 2 𝛿𝛿𝐴𝐴 [𝑋𝑋,𝜉𝜉 𝑛𝑛 ]⎟ 𝛿𝛿𝛾𝛾21 ⎟ 0 ⎟ 0 ⎟ 2 𝛿𝛿𝐴𝐴 [𝑋𝑋,𝜉𝜉 𝑛𝑛 ]⎟ 𝛿𝛿𝑉𝑉23 ⎟ 𝛿𝛿𝐴𝐴2 [𝑋𝑋,𝜉𝜉 𝑛𝑛 ]⎟ 𝛿𝛿𝑉𝑉24 ⎟ 0 ⎟ 𝛿𝛿𝛼𝛼𝑜𝑜𝑜𝑜𝑜𝑜 𝛿𝛿𝐴𝐴2 [𝑋𝑋,𝜉𝜉 𝑛𝑛 ] 𝛿𝛿𝑉𝑉21 ⎠ Оператор (2.2.7) представляет собой матрицу частных производных модели по параметрам, где количество столбцов определяется количеством элементов в потоке входных данных X, а количество строк равно количеству параметров в векторе 𝜉𝜉 𝑛𝑛 . Процесс заканчивается, когда разница между реальным и модельным дебитами 𝜑𝜑 𝑛𝑛 принимает удовлетворительно малое значение, либо не удалось получить новый вектор 𝜉𝜉 𝑛𝑛 Проверка на тестовом примере. В качестве примера для проверки адекватности поведения модели были использованы данные по эксплуатации тестового месторождения, состоящего из 4 добывающих скважин (№ 1, 2, 3, 4) и 5 нагнетательных скважин (№ 5, 6, 7, 8, 9). 40

Сетка скважин месторождения По результатам решения оптимизационной задачи для экспериментального месторождения модельные дебиты приблизились к дебитам из истории разработки с относительной погрешностью 7,8%. Временной ряд истории разработки тестового месторождения составил 50 месяцев, из которых 40 было выделено на настройку модели, поэтому период прогноза составил 10 оставшихся месяцев. На основе рассчитанных модельных параметров был сделан прогноз динамики добычи жидкости по скважинам на определенный период прогноз. На графиках представлены результаты прогноза по каждой добывающей скважине (Рисунок 2.2.2). а) скважина № 1 б) скважина № 2 41

в) скважина № 3 г) скважина № 4 Сравнение модельного (пунктирная линия) дебита и дебита из истории разработки тестового месторождения Относительная погрешность прогноза по каждой скважине составила 7,42%, 8,65%, 6,04% и 9,1% соответственно. Такие результаты прогноза показывают, что построенная модель ведет себя адекватно и ее можно использовать для томографических вычислений путем моделирования депрессии в нагнетательных скважинах и регистрации отклика на эту депрессию в добывающих скважинах. На графиках (Рисунок 2.2.3) представлены результаты измерения времени движения сигнала между одной из нагнетательных скважин (№5) и всеми добывающими скважинами (№1, 2, 3, 4). Время прохождения сигнала между нагнетательной скважиной и добывающими скважинами 42

В дальнейшем полученные интервальные времена, будут использоваться в качестве начальных данных в задаче гидродинамической томографии для реконструкции пространственного распределения фильтрационного сопротивления проницаемого пласта нефтегазовых месторождений, технология и вычислительная схема которой представлена в 3 главе. 2.3. Заключение Разработанная модель может быть применена для расчета синтетических временных интервалов движения экстремальной точки между всеми парами скважин адекватен модели динамики эксплуатации залежи и обеспечивает получение входных данных для технологии гидродинамической томографии. Эти данные могут быть использованы для реконструкции пространственного распределения фильтрационного сопротивления в пределах гидродинамически связанных компонент проницаемого пласта. Итерационный процесс подбора оптимального вектора параметров для математической модели динамики эксплуатации нефтегазового месторождения соответствует первому защищаемому положению: − математической метод моделирования и вычислительная схема пассивной гидродинамической томографии для синтезирования данных с целью включения их в активную форму томографии. 43

ГЛАВА 3. МЕТОДИКА И ТЕХНОЛОГИЯ РЕШЕНИЯ ЗАДАЧИ ГИДРОДИНАМИЧЕСКОЙ ТОМОГРАФИИ ПРОНИЦАЕМОГО ПЛАСТА В представленной главе диссертационной работы описаны принципы и методы движения экстремальной точки (точка перегиба на КВД), разработаны алгоритмы и вычислительные схемы построения прогноза параметров проницаемости нефтегазового пьезопроводности по пласта, данным выраженные через гидродинамического коэффициент прослушивания, преобразованного в томографическую форму за счет системной организации данных. Вычислительная основа методов гидродинамической томографии адаптирована к уравнениям движения потока флюидов в неоднородной среде. 3.1. Организация и обработка результатов гидродинамического прослушивания скважин в томографическую систему данных Проблема поиска и определения пространственного распределения фильтрационного сопротивления движения потоков флюидов внутри продуктивных пластов нефтегазовых месторождений была представлена как томографическая задача на начальном этапе в работах [33, 45, 101] и зарегистрированном патенте «Способ гидродинамической томографии проницаемости пласта» [110]. Метод гидродинамической томографии состоит в построении томографической системы наблюдений за гидродинамическим прослушиванием скважин и обработке полученных веерных измерений интервальных времен прихода реакции в виде изменения уровня давления в скважинах-приемниках, на создаваемую депрессию в скважинах-источниках. Пусть 𝑆𝑆 – проекция на поверхность исследуемого продуктивного пласта нефтегазового месторождения. На проекцию 𝑆𝑆 накладывается регулярная сетка, в ячейках которой определены 𝑁𝑁 скважин, которые используются последовательно в качестве источника и приемника создаваемой депрессии на пласт с координатами на сети 𝜉𝜉𝑖𝑖 = {𝑥𝑥𝑖𝑖 , 𝑦𝑦𝑖𝑖 } ∈ 𝑆𝑆, 𝑖𝑖 = 1 ÷ 𝑁𝑁. 44

Из фонда всех (или части) добывающих и нагнетательных скважин последовательно определяется скважина 𝜉𝜉𝑖𝑖 , которая выбирается как источник депрессии и в которую производится закачка жидкости с некоторым постоянным уровнем дебита. Далее согласно принципах реализации гидродинамического прослушивания скважин изменяется внутрипластовое давление и на окружающих гидродинамически связанных скважинах 𝜉𝜉𝑗𝑗 с источником (депрессии) регистрируется наблюдается время наступления изменения отклика на режима созданную работы и депрессию, связанный с действием распространения фронта давления в проницаемом пласте нефтегазового месторождения. Регистрируются интервалы времен с момента начала депрессии на скважине-источнике и пике изменения давления в скважине-приемнике в форме кривой восстановления (изменения) давления (КВД). На кривой восстановления давления определяется особая (некоторая реперная) точка. Рассчитывается интервал времени 𝜏𝜏(𝜉𝜉𝑖𝑖 , 𝜉𝜉𝑗𝑗 ) ее перемещение между парой скважин от скважины с номером 𝑖𝑖 к скважине с номером 𝑗𝑗. В качестве такой точки Щелкачевым В. Н. впервые было предложено рассматривать точку перегиба на КВД [23, 27]. Объединение данных гидродинамического прослушивания скважин нефтегазового месторождения в томографическую систему обработки данных реализуется через следующий алгоритм. В системе 𝑁𝑁 скважин определяется скважина-источник депрессии (нагнетательная) с номером 𝑖𝑖, находящаяся в границе продуктивного пласта гидродинамически связанная с другими скважинами и подвергается длительной депрессии (депрессия на пласт возможна как положительная, так и отрицательная). В скважинах-приемниках депрессии 𝑀𝑀𝑖𝑖 < 𝑁𝑁 (всех или части доступных для эксперимента скважин) определяется интервал времени 𝜏𝜏𝑖𝑖𝑖𝑖 , 𝑗𝑗 = 1 ÷ 𝑀𝑀𝑖𝑖 , в течение которого депрессия достигла скважин 𝑀𝑀𝑖𝑖 и установила максимальную скорость отдачи на скважине, т. е. была достигнута точки перегиба нарастания (или спада) давления на КВД. 45

Всего 𝑀𝑀𝑖𝑖 значений интервальных времен для веера распространения фронта давления в виде точек перегиба на кривых изменения давления от скважиныисточника 𝑖𝑖 к доступным скважинам-приемникам. Далее выбирается новая скважина-источник депрессии 𝑖𝑖 + 1, процесс повторяется и рассчитывается новый следующий веер значений 𝜏𝜏𝑖𝑖+1,𝑗𝑗 , 𝑗𝑗 = 1 ÷ 𝑀𝑀𝑖𝑖+1 . Эксперимент продолжается для максимально возможного количества скважин, в итоге получаем матрицу значений интервальных времен движения реперной точки между всеми парами скважин по площади проницаемого пласта нефтегазового месторождения. В этом наборе полученных измерений интервальных времен наблюдается типичная томографическая система, путем наложения результатов нескольких отдельных экспериментов друг на друга для создания единой карты [15, 50]. По результатам предложенных экспериментальных измерений реализуется задача прогнозирования (или реконструкции) пространственного распределения фильтрационного коэффициент пьезопроводности, сопротивления, который выраженного характеризует через процесс распространения давления в проницаемой системе. Реконструкция базируется на способах решения некорректных задач для уравнений интервалов времени движения особых точек восстановления давлений, между выделенными парами скважин. Для реализации предложенной реконструкции необходимо: построить интеграл, характеризующий уравнения полученной системы интервальных времен; разработать приближенный способ моделирования для построенного интеграла; для найденных уравнений разработать примерные устойчивые методы решения обратных задач. 3.2. Алгоритм расчета интервальных времен движения особых точек динамики восстановления давления В работе [66] представлено несколько примеров задач по нахождению оптимального пути между двумя точками двумерного пространства, выходящей из точки 𝑥𝑥0 и заканчивающийся в точке 𝑥𝑥. Траектория движения 46

описывается некоторой кривой 𝑥𝑥(𝑠𝑠). Существует критерий 𝐽𝐽(𝑥𝑥(𝑠𝑠)), который оценивает оптимальность рассматриваемого участка траектории. Для нахождения оптимальной траектории решается задача: 𝐽𝐽 �𝑥𝑥(𝑠𝑠 𝑒𝑒𝑒𝑒𝑒𝑒 )� → 𝑚𝑚𝑚𝑚𝑚𝑚 � 𝑥𝑥(0) = 𝑥𝑥0 В работе рассматривается пример движения особой точки вдоль линии луча 𝑥𝑥(𝑠𝑠) в поле координат, подчиняющейся принципу Ферма, который утверждает, что время перемещения точки между началом и окончанием пути должно определяться как минимальное, т. е. давление в среде распространяется по пути наименьшего сопротивления проницаемой среды [30]: 𝑥𝑥 𝑑𝑑𝑑𝑑 → 𝑚𝑚𝑖𝑖𝑛𝑛 𝑥𝑥0 𝑉𝑉(𝑥𝑥(𝑠𝑠)) 𝜏𝜏 = � Решение поставленной задачи (3.2.2) реализуется по вышеобозначему принципу, где оптимизируется (минимизируется) интеграл движения особой точки перегиба на рассчитывается кривой на основе восстановления давления принципа Беллмана во времени и динамического программирования: 𝑙𝑙(𝜉𝜉)𝑑𝑑𝑑𝑑 � = 𝜏𝜏(𝑥𝑥0 , 𝑥𝑥) 3𝜅𝜅(𝜉𝜉) 𝐿𝐿(𝑥𝑥0 ,𝑥𝑥) min �� 𝐿𝐿(𝑥𝑥0 ,𝑥𝑥) Начальная 𝑥𝑥0 и конечная 𝑥𝑥 точки траектории движения определены, как координаты скважин для которых рассчитывается оптимальный интервал времени движения. Характерной особенностью(3.2.3) от (3.2.2) является вхождение оптимальной траектории пути в определяемый интеграл как траектории интегрирования. То есть интегрируемое выражение является функциональной зависимостью «скорости движения луча» от длинны пройденной траектории между точками 𝑉𝑉(𝜉𝜉) = 3𝜅𝜅(𝜉𝜉) 𝑙𝑙(𝜉𝜉) . Из-за того, что скорость движения уменьшается пропорционально длине пройденной траектории, это вызывает существенные изменения кинематики перемещения особой точки в проницаемой среде. Но все же, основные 47 принципы использования

Беллмановского алгоритма при решении задачи динамического программирования также распространяются и на (3.2.2). В данном разделе представлена реализация решения прямой задачи гидродинамической томографии проницаемого пласта на основе поиска оптимальной (кратчайшей) траектории пути между парами скважин для исследуемой базирующаяся математической на принципе модели нефтегазового динамического месторождения, программирования [58]. С помощью колотого выполняется решение некорректной задачи, основанное на итерационном подходе последовательного уточнения решения [30]. Выбрана некоторая область 𝑆𝑆 в пределах исследуемого проницаемого пласта нефтегазового месторождения, на которой расположена сетка скважин и выполняются необходимые построения. Координата точки на области обозначается 𝜉𝜉, которая покрыта регулярной сеткой так, что каждая точка области в сеточном представлении однозначно определяется парой индексов 𝜉𝜉(𝑖𝑖, 𝑗𝑗) = (𝑖𝑖, 𝑗𝑗). На сети нулевое распределение коэффициента пьезопроводности 𝜅𝜅(𝜉𝜉) = 𝜅𝜅(𝑖𝑖, 𝑗𝑗) может быть задано произвольным образом или на основе имеющейся информации о характере распределения фильтрационно-емкостных свойств в пределах исследуемого месторождения (см. рисунок 3.2.1) участка. На сетке скважин рассматриваются пары скважин 𝑞𝑞(𝑚𝑚, 𝑛𝑛), 𝑚𝑚 = 1 ÷ 𝑀𝑀, 𝑛𝑛 = 1 ÷ 𝑁𝑁, где 𝑚𝑚 – номер скважины-источника депрессии на пласт и 𝑛𝑛 – номер скважины-приемника возмущения. Всего в эксперименте принимает 𝑄𝑄 пар скважин: 𝑄𝑄 = 𝑀𝑀 × 𝑁𝑁. 48

Исходные данные для расчета оптимального пути Распространение депрессии в среде между парой скважин совершается по линиям траектории пути 𝐿𝐿𝑞𝑞 = 𝐿𝐿(𝜉𝜉𝑞𝑞𝑞𝑞 , 𝜉𝜉𝑞𝑞𝑞𝑞 ), где начало траектории в точке 𝜉𝜉𝑞𝑞𝑞𝑞 = (𝑖𝑖, 𝑗𝑗)𝑞𝑞𝑞𝑞 и заканчивающейся в точке 𝜉𝜉𝑞𝑞𝑞𝑞 = (𝑖𝑖, 𝑗𝑗)𝑞𝑞𝑞𝑞 . Траектории пронумерованы индексами 𝑞𝑞 = 1 ÷ 𝑄𝑄, интервальное время движения возмущения по траектории 𝐿𝐿𝑞𝑞 есть 𝜏𝜏𝑞𝑞 . Количество пройденных точек сетки в каждой траектории есть 𝑞𝑞𝑞𝑞 = 1 ÷ 𝑁𝑁𝑞𝑞 . Текущая координата на траектории 𝐿𝐿𝑞𝑞 определяется переменной 𝜉𝜉𝑞𝑞 = (𝑖𝑖, 𝑗𝑗)𝑞𝑞 = (𝑖𝑖, 𝑗𝑗, 𝑞𝑞𝑞𝑞)𝑞𝑞 . Интервальное время распространения движения волны 𝜏𝜏𝑞𝑞 по 𝐿𝐿𝑞𝑞 определяется по формуле: 𝑁𝑁𝑞𝑞 � 𝑞𝑞𝑞𝑞=1 𝑙𝑙(𝑞𝑞𝑞𝑞) = 𝜏𝜏𝑞𝑞 , 3𝜅𝜅(𝑞𝑞𝑞𝑞) где 𝑙𝑙(𝑞𝑞𝑞𝑞) – длина части траектории 𝐿𝐿𝑞𝑞 до ячейки 𝑞𝑞𝑞𝑞. Наличие этой компоненты проявляется в уменьшении скорости движения по мере продвижения особой точки вдоль траектории пути: 𝑉𝑉 = 3𝜅𝜅(𝑞𝑞𝑞𝑞) 𝑙𝑙(𝑞𝑞𝑞𝑞) 49

Траектория 𝐿𝐿𝑞𝑞 , соединяющая точки 𝜉𝜉𝑚𝑚 и 𝜉𝜉𝑛𝑛 есть линия, служащая решением задачи: min � 𝐺𝐺(𝜉𝜉𝑚𝑚 ,𝜉𝜉𝑛𝑛 ) � 𝐺𝐺(𝜉𝜉𝑚𝑚 ,𝜉𝜉𝑛𝑛 ) 𝑙𝑙(𝜉𝜉) 𝑑𝑑𝜉𝜉 � = 3𝜅𝜅(𝜉𝜉) = 𝜏𝜏𝑞𝑞 , 𝑁𝑁𝑞𝑞 � 𝐿𝐿𝑞𝑞 (𝜉𝜉𝑚𝑚 ,𝜉𝜉𝑛𝑛 ) 𝑙𝑙(𝜉𝜉) 𝑙𝑙(𝑞𝑞𝑞𝑞) 𝑑𝑑𝜉𝜉 ≈ � 3𝜅𝜅(𝜉𝜉) 3𝜅𝜅(𝑞𝑞𝑞𝑞) 𝑞𝑞𝑞𝑞=1 где 𝐺𝐺(𝜉𝜉𝑚𝑚 , 𝜉𝜉𝑛𝑛 ) – возможные траектории между парами скважин 𝑞𝑞(𝑚𝑚, 𝑛𝑛). В итоге более длинный путь, но который обходит участки повышенных значений 𝜅𝜅(𝑞𝑞𝑞𝑞) коэффициента пьезопроводности, за счет уменьшения скорости движения по гиперболическому закону может оказаться проигрышным в сравнении с более коротким, но который пересекает локальные максимумы 𝜅𝜅(𝑞𝑞𝑞𝑞) на карте. Это специфика динамики распространения особой точки коэффициента пьезопроводности. Способ вычисления траектории 𝐿𝐿𝑞𝑞 (𝑚𝑚, 𝑛𝑛), по которой определяются значения 𝑙𝑙(𝑞𝑞𝑘𝑘), основан на последовательном нахождении окончания траектории начиная с текущего индекса 𝑧𝑧𝑧𝑧𝑞𝑞 для траектории 𝑞𝑞 до 𝑞𝑞𝑞𝑞 вместе с «остатком траектории» 𝐿𝐿𝑞𝑞 (𝑞𝑞𝑞𝑞), 𝑞𝑞𝑞𝑞 = 𝑧𝑧𝑧𝑧𝑞𝑞 ÷ 𝑞𝑞𝑞𝑞 в предположении, что 𝐿𝐿𝑞𝑞 (𝑞𝑞𝑞𝑞), 𝑞𝑞𝑞𝑞 = 1 ÷ 𝑧𝑧𝑧𝑧𝑞𝑞 и, следовательно, 𝑙𝑙(𝑞𝑞𝑞𝑞) заданы своими нулевыми приближениями: 𝐿𝐿𝑞𝑞 0 (𝑞𝑞𝑞𝑞), 𝑙𝑙𝑞𝑞 0 (𝑞𝑞𝑞𝑞), 𝑞𝑞𝑞𝑞 = 1 ÷ 𝑧𝑧𝑧𝑧𝑞𝑞 . min 𝑧𝑧𝑘𝑘𝑞𝑞 𝑞𝑞𝑞𝑞 𝑙𝑙𝑞𝑞 0 (𝑞𝑞𝑞𝑞) 𝑔𝑔𝑞𝑞 (𝑞𝑞𝑞𝑞) + � �� = 𝜏𝜏𝑞𝑞 3𝜅𝜅(𝑞𝑞𝑞𝑞) 3𝜅𝜅(𝑞𝑞𝑞𝑞) 𝑔𝑔(𝑞𝑞𝑞𝑞), 𝑞𝑞𝑞𝑞=𝑧𝑧𝑘𝑘𝑞𝑞 ÷𝑞𝑞𝑞𝑞 𝑞𝑞𝑞𝑞=1 𝑞𝑞𝑞𝑞=𝑧𝑧𝑘𝑘𝑞𝑞 где 𝑔𝑔𝑞𝑞 (𝑞𝑞𝑞𝑞) – подобранные длины второй половины траектории для 𝑞𝑞𝑞𝑞 = 𝑧𝑧𝑘𝑘𝑞𝑞 ÷ 𝑞𝑞𝑞𝑞. Алгоритм расчета реализован на последовательном поиске решения и сравнении результатов решения задачи, начиная от 𝑧𝑧𝑘𝑘𝑞𝑞 = 𝑞𝑞𝑞𝑞 − 1 до 𝑧𝑧𝑘𝑘𝑞𝑞 = 𝑞𝑞𝑞𝑞 + 1. При этом для каждого этапа 𝑞𝑞𝑞𝑞 сохраняются для последующего сравнения все значения конца траектории 𝑔𝑔𝑞𝑞 (𝑧𝑧𝑘𝑘𝑞𝑞 ), служащие именем завершающего участка 𝑔𝑔𝑞𝑞 (𝑞𝑞𝑞𝑞), 𝑞𝑞𝑞𝑞 = 𝑧𝑧𝑧𝑧𝑞𝑞 ÷ 𝑞𝑞𝑞𝑞. 50

После завершения счета, полученная траектория движения особой точки 𝐺𝐺𝑞𝑞 (𝑞𝑞𝑞𝑞), характеризующаяся длинами 𝑔𝑔𝑞𝑞 (𝑞𝑞𝑞𝑞), принимается за 𝐿𝐿1 (𝑞𝑞𝑞𝑞) и 𝑙𝑙1 (𝑞𝑞𝑞𝑞) соответственно. Рассмотрим скважины с координатами (𝑖𝑖𝑚𝑚 , 𝑗𝑗𝑚𝑚 ) и (𝑖𝑖𝑛𝑛 , 𝑗𝑗𝑛𝑛 ). Согласно принципу Беллмана, движение при расчете траекторий начинается с координаты скважина-источника, т.е. с точки 𝜉𝜉𝑛𝑛 = 𝜉𝜉(𝑖𝑖𝑛𝑛 , 𝑗𝑗𝑛𝑛 ) и заканчивается в 𝜉𝜉𝑚𝑚 = 𝜉𝜉(𝑖𝑖𝑚𝑚 , 𝑗𝑗𝑚𝑚 ). Этап 1 (рисунок 3.2.2). На первом этапе 𝑧𝑧𝑘𝑘𝑞𝑞 = 𝑞𝑞𝑞𝑞 − 1 ищем 𝐽𝐽 путей 𝐺𝐺𝑞𝑞 1 (𝑞𝑞𝑞𝑞) из 𝜉𝜉𝑛𝑛 = 𝜉𝜉(𝑖𝑖𝑛𝑛 , 𝑗𝑗𝑛𝑛 ) в ячейки (𝑖𝑖𝑛𝑛 − 1, 𝑗𝑗1 ), 𝑗𝑗1 = 1 ÷ 𝐽𝐽 и столько же путей 𝐿𝐿𝑞𝑞 0 (𝑞𝑞𝑞𝑞), 𝑙𝑙𝑞𝑞 0 (𝑞𝑞𝑞𝑞), 𝑞𝑞𝑞𝑞 = 1 ÷ 𝑧𝑧𝑧𝑧𝑞𝑞 из (𝑖𝑖𝑚𝑚 , 𝑗𝑗𝑚𝑚 ) в (𝑖𝑖𝑛𝑛 − 1, 𝑗𝑗1 ). На первом этапе получаем следующий путь: 𝐿𝐿𝑞𝑞 1 (𝑚𝑚, 𝑛𝑛) = 𝐿𝐿𝑞𝑞 0 (𝑞𝑞𝑞𝑞) + 𝐺𝐺𝑞𝑞 1 (𝑞𝑞𝑞𝑞) и соответствующие им длины: 𝑙𝑙𝑞𝑞 1 (𝑚𝑚, 𝑛𝑛) = 𝑙𝑙𝑞𝑞 0 (𝑞𝑞𝑞𝑞) + 𝑔𝑔𝑞𝑞 1 (𝑞𝑞𝑞𝑞) Этап 2 (рисунок 3.2.3). На этом этапе 𝑧𝑧𝑘𝑘𝑞𝑞 = 𝑞𝑞𝑞𝑞 − 2. Далее из 𝜉𝜉𝑛𝑛 = 𝜉𝜉(𝑖𝑖𝑛𝑛 , 𝑗𝑗𝑛𝑛 ) ищем пути 𝐺𝐺𝑞𝑞 2 (𝑞𝑞𝑞𝑞) в ячейки (𝑖𝑖𝑛𝑛 − 2, 𝑗𝑗2 ), через ячейки, для которых были найдены пути на этапе 1 (𝑖𝑖𝑛𝑛 − 1, 𝑗𝑗1 ), 𝑗𝑗2 = 1 ÷ 𝐽𝐽. И столько же путей из (𝑖𝑖𝑚𝑚 , 𝑗𝑗𝑚𝑚 ) в (𝑖𝑖𝑛𝑛 − 2, 𝑗𝑗2 ). На рисунке 3.2.3 показан путь для одной из 𝐽𝐽 ячеек. 51

Графическое представление этапа 1 Графическое представление этапа 2 Этап 3 (рисунок 3.2.4). Рассчитываем оптимальные траектории 𝐿𝐿𝑞𝑞 2 (𝑚𝑚, 𝑛𝑛) = 𝐿𝐿𝑞𝑞 0 (𝑞𝑞𝑞𝑞) + 𝐺𝐺𝑞𝑞 2 (𝑞𝑞𝑞𝑞) и длины 𝑙𝑙𝑞𝑞 2 (𝑚𝑚, 𝑛𝑛) по условию, в формуле (3.2.7). На рисунке 52

3.2.4 показан найденный путь, проходящий через ячейку (𝑖𝑖𝑛𝑛 − 2, 𝑗𝑗2 ), из ячейки 𝜉𝜉𝑛𝑛 = 𝜉𝜉(𝑖𝑖𝑛𝑛 , 𝑗𝑗𝑛𝑛 ) в ячейку (𝑖𝑖𝑚𝑚 , 𝑗𝑗𝑚𝑚 ). Пути 𝐺𝐺𝑞𝑞 2 (𝑞𝑞𝑞𝑞)𝑜𝑜𝑜𝑜𝑜𝑜 запоминаем, далее работаем с ними. Этап 4. Если 𝑖𝑖𝑛𝑛 − 2 ≠ 𝑖𝑖𝑚𝑚 и 𝑗𝑗2 ≠ 𝑗𝑗𝑚𝑚 , перейти к шагу 5. В противном случае перейти к шагу 8. Этап 5 (рисунок 3.2.5). Если мы уже на итерации 𝑧𝑧𝑘𝑘𝑞𝑞 = 𝑞𝑞𝑞𝑞 − 𝑧𝑧. Далее из 𝜉𝜉𝑛𝑛 = 𝜉𝜉(𝑖𝑖𝑛𝑛 , 𝑗𝑗𝑛𝑛 ) находим пути 𝐺𝐺𝑞𝑞 𝑧𝑧 (𝑞𝑞𝑞𝑞) в ячейки (𝑖𝑖𝑛𝑛 − 𝑧𝑧, 𝑗𝑗𝑧𝑧 ), через ячейки, в которых были найдены пути на 𝑧𝑧 − 1 шаге (𝑖𝑖𝑛𝑛 − 𝑧𝑧 − 1, 𝑗𝑗𝑧𝑧−1 ), 𝑗𝑗𝑧𝑧 = 1 ÷ 𝐽𝐽. И столько же путей из (𝑖𝑖𝑚𝑚 , 𝑗𝑗𝑚𝑚 ) в (𝑖𝑖𝑛𝑛 − 𝑧𝑧, 𝑗𝑗𝑧𝑧 ). Этап 6. C учетом формулы (3.2.6), как условия оптимальности пути, находим 𝐿𝐿𝑞𝑞 𝑁𝑁𝑞𝑞 (𝑚𝑚, 𝑛𝑛) = 𝐿𝐿𝑞𝑞 0 (𝑞𝑞𝑞𝑞) + 𝐺𝐺𝑞𝑞 𝑧𝑧 (𝑞𝑞𝑞𝑞)𝑜𝑜𝑜𝑜𝑜𝑜 и длины 𝑙𝑙𝑞𝑞 𝑧𝑧 (𝑚𝑚, 𝑛𝑛). На рисунке 3.2.7 показан путь для одной из 𝐽𝐽 ячеек при 𝑧𝑧 = 3. Графическое представление этап 3 53

Графическое представление этап 5 для 𝑧𝑧 = 3 Шаг 8 54

Выбор оптимальной траектории движения точки кривой восстановления давления Этап 7. Если 𝑖𝑖𝑛𝑛 − 𝑧𝑧 ≠ 𝑖𝑖𝑚𝑚 и 𝑗𝑗𝑧𝑧 ≠ 𝑗𝑗𝑚𝑚 , перейти к этапу 5. В противном случае перейти к этапу 8. Этап 8 (рисунки 3.2.6, 3.2.7). На последней итерации 𝑧𝑧𝑘𝑘𝑞𝑞 = 𝑞𝑞𝑞𝑞 + 1 рассчитываем 𝐽𝐽 значений 𝐿𝐿𝑞𝑞 (𝑚𝑚, 𝑛𝑛) = 𝐿𝐿𝑞𝑞 0 (𝑞𝑞𝑞𝑞) + 𝐺𝐺𝑞𝑞 𝑧𝑧 (𝑞𝑞𝑞𝑞). Далее находим минимальное по 𝑗𝑗𝑧𝑧 значение интервального времени распространения возмущения между скважинами, расположенными в точках 𝜉𝜉𝑛𝑛 = 𝜉𝜉(𝑖𝑖𝑛𝑛 , 𝑗𝑗𝑛𝑛 ) и 𝜉𝜉𝑚𝑚 = 𝜉𝜉(𝑖𝑖𝑚𝑚 , 𝑗𝑗𝑚𝑚 ) соответствующее значению индекса 𝑗𝑗𝑧𝑧 . Оптимальная траектория, удовлетворяющая условию (3.3.7) и соединяющая точки (𝑖𝑖𝑛𝑛 , 𝑗𝑗𝑛𝑛 ) и (𝑖𝑖𝑚𝑚 , 𝑗𝑗𝑚𝑚 ) есть: 3.3. 𝐿𝐿𝑞𝑞 (𝑚𝑚, 𝑛𝑛) = min[𝐿𝐿𝑞𝑞 0 (𝑞𝑞𝑞𝑞) + 𝐺𝐺𝑞𝑞 𝑧𝑧 (𝑞𝑞𝑞𝑞)𝑜𝑜𝑜𝑜𝑜𝑜 ] 𝑗𝑗𝑧𝑧 Вычислительная схема метода гидродинамической томографии Рассмотрим поэтапно вычислительную схему гидродинамической томографии [48], после выполнения всех действий по расчету оптимальных 55

интервальных времен 𝜏𝜏𝑞𝑞 распространения движения волны внутри пласта, соответствующих заданному пространственному распределению фильтрационного сопротивления. Время распространения движения волны 𝜏𝜏𝑞𝑞 рассчитывается при помощи следующего оператора: 𝐴𝐴[𝜅𝜅(𝜉𝜉)] = � 𝐿𝐿𝑞𝑞 𝑙𝑙(𝜉𝜉)𝑑𝑑𝑑𝑑 = 𝜏𝜏𝑞𝑞 3𝜅𝜅(𝜉𝜉) что в дискретном представлении: (𝑖𝑖,𝑗𝑗)=(𝑖𝑖,𝑗𝑗)𝑞𝑞𝑞𝑞 (𝑖𝑖,𝑗𝑗)=(𝑖𝑖,𝑗𝑗)𝑞𝑞𝑞𝑞 𝑙𝑙(𝑖𝑖, 𝑗𝑗) = 𝜏𝜏𝑞𝑞 3𝜅𝜅(𝑖𝑖, 𝑗𝑗) восстановления давления), 𝐴𝐴̅[𝜅𝜅(𝜉𝜉)] = � Пусть 𝝉𝝉� = {𝜏𝜏̅𝑞𝑞 } есть наблюдаемое время распространения давления (точки перегиба кривой эксперимента проведенного на измеренное математической в результате модели пассивной гидродинамической томографии и которое соответствует неизвестному пространственному распределению коэффициента пьезопроводности 𝜅𝜅(𝜉𝜉) + ∆𝜅𝜅(𝜉𝜉). Задача состоит в том, чтобы определить ∆𝜅𝜅(𝜉𝜉) такое, что рассчитанное 𝜏𝜏𝑞𝑞 было максимально приближено к 𝜏𝜏̅𝑞𝑞 . Тогда разность ∆𝜏𝜏𝑞𝑞 наблюдаемого и рассчитанного интервальных времен получаем в виде соотношения, которое выражается из нескольких приближенных равенств: 𝜏𝜏̅𝑞𝑞 − 𝜏𝜏𝑞𝑞 = � 𝐿𝐿𝑞𝑞 � � 𝐿𝐿𝑞𝑞 𝑙𝑙(𝜉𝜉) 𝑙𝑙(𝜉𝜉) 𝑑𝑑𝑑𝑑 − � 𝑑𝑑𝑑𝑑 ≈ 3[𝜅𝜅(𝜉𝜉) + ∆𝜅𝜅(𝜉𝜉)] 𝐿𝐿𝑞𝑞 3𝜅𝜅(𝜉𝜉) 𝑙𝑙(𝜉𝜉) 𝑙𝑙(𝜉𝜉) − � 𝑑𝑑𝑑𝑑 3[𝜅𝜅(𝜉𝜉) + ∆𝜅𝜅(𝜉𝜉)] 3𝜅𝜅(𝜉𝜉) Меняя знак, получаем: A′ [𝜅𝜅(𝜉𝜉)]∆𝜅𝜅(𝜉𝜉) = 𝜏𝜏𝑞𝑞 − 𝜏𝜏̅𝑞𝑞 = � 𝐿𝐿𝑞𝑞 𝑙𝑙(𝜉𝜉)∆𝜅𝜅(𝜉𝜉) 𝑑𝑑𝑑𝑑 = ∆𝜏𝜏𝑞𝑞 𝑑𝑑𝑑𝑑 3𝜅𝜅 2 (𝜉𝜉) Оператор A′ [𝜅𝜅(𝜉𝜉)] суть производная Фреше к оператору 𝐴𝐴[𝜅𝜅(𝜉𝜉)] в «точке» 𝜅𝜅(𝜉𝜉). A′ [𝜅𝜅(𝜉𝜉)] – линейный оператор, действующий на вектор ∆𝜅𝜅(𝜉𝜉). Уравнения (3.3.3) и (3.3.4) получены в предположении того, что величина 56

∆𝜅𝜅(𝜉𝜉) настолько мала, что, произведение (𝜅𝜅(𝜉𝜉) + ∆𝜅𝜅(𝜉𝜉)) × 𝜅𝜅(𝜉𝜉) ≈ 𝜅𝜅(𝜉𝜉)2 и траектории 𝐿𝐿𝑞𝑞 = 𝐿𝐿(𝜉𝜉𝑞𝑞𝑞𝑞 , 𝜉𝜉𝑞𝑞𝑞𝑞 ) для пьезопроводности 𝜅𝜅(𝜉𝜉) + ∆𝜅𝜅(𝜉𝜉) и 𝜅𝜅(𝜉𝜉) с достаточной точностью совпадают. Уравнение (3.3.4) в дискретной форме выглядит следующим образом: �′ [𝜅𝜅(𝜉𝜉)]∆𝜅𝜅(𝜉𝜉) = A (𝑖𝑖,𝑗𝑗)=(𝑖𝑖,𝑗𝑗)𝑞𝑞𝑞𝑞 � (𝑖𝑖,𝑗𝑗)=(𝑖𝑖,𝑗𝑗)𝑞𝑞𝑞𝑞 𝑙𝑙(𝑖𝑖, 𝑗𝑗)∆𝜅𝜅(𝑖𝑖, 𝑗𝑗) = ∆𝜏𝜏𝑞𝑞 3𝜅𝜅 2 (𝑖𝑖, 𝑗𝑗) Если количество узлов сети в траектории 𝑞𝑞 равно 𝑁𝑁𝑞𝑞 , то размерность 𝑄𝑄 вектора ∆𝜅𝜅(𝑖𝑖, 𝑗𝑗) равна С = ∑𝑞𝑞=1 𝑁𝑁𝑞𝑞 . Областью значения преобразования (3.3.5) служит вектор ∆𝝉𝝉 = {∆𝜏𝜏𝑞𝑞 } размерности 𝑄𝑄. После того как вектор ∆𝜅𝜅(𝑖𝑖, 𝑗𝑗) найден, следующее приближение к распределению коэффициента пьезопроводности определяется правилом: 𝜅𝜅1 (𝜉𝜉) = 𝜅𝜅(𝜉𝜉) + ∆𝜅𝜅(𝜉𝜉) Вычислительной основой [49] для решения задачи (3.3.5) относительно вектора ∆𝜅𝜅(𝑖𝑖, 𝑗𝑗) в точках 𝜉𝜉(𝑖𝑖, 𝑗𝑗) траекторий 𝐿𝐿𝑞𝑞 = 𝐿𝐿(𝜉𝜉𝑞𝑞𝑞𝑞 , 𝜉𝜉𝑞𝑞𝑞𝑞 ) служит итерационный процесс, на 𝑧𝑧 + 1 итерации имеющий вид: ∆𝜅𝜅 𝑧𝑧+1 (𝑖𝑖, 𝑗𝑗) = ∆𝜅𝜅 𝑧𝑧 (𝑖𝑖, 𝑗𝑗) + 𝛼𝛼 𝑧𝑧 𝐴𝐴′∗ [𝜅𝜅(𝜉𝜉)]𝝋𝝋𝒛𝒛 � ||∆𝜅𝜅 𝑧𝑧+1 − ∆𝜅𝜅 𝑧𝑧 || ≤ 𝜀𝜀 𝝋𝝋𝒛𝒛 = 𝐴𝐴[𝜅𝜅(𝜉𝜉) + ∆𝜅𝜅 𝑧𝑧 (𝑖𝑖, 𝑗𝑗)] − 𝝉𝝉� = �� 𝐿𝐿𝑞𝑞 𝑙𝑙(𝜉𝜉) 𝑑𝑑𝑑𝑑 − 𝜏𝜏̅𝑞𝑞 , 𝑞𝑞 = 1 ÷ 𝑄𝑄� 3[𝜅𝜅(𝜉𝜉) + ∆𝜅𝜅 𝑧𝑧 (𝑖𝑖, 𝑗𝑗)] где 𝝋𝝋𝒛𝒛 – разница между наблюдаемыми интервальными временами и временами, рассчитанными на итерации 𝑧𝑧. Нулевое приближение принимаем равным ∆𝜅𝜅 0 (𝑖𝑖, 𝑗𝑗) = 0. Числовой параметр 𝛼𝛼 𝑧𝑧 – параметр релаксации, подбираемый на каждом шаге так, чтобы итерационный процесс (3.3.7) сходился: ‖𝐴𝐴′∗ [𝜅𝜅(𝜉𝜉)]𝝋𝝋𝒛𝒛 ‖2 𝛼𝛼 = ′ ‖𝐴𝐴 [𝜅𝜅(𝜉𝜉)]𝐴𝐴′∗ [𝜅𝜅(𝜉𝜉)]𝝋𝝋𝒛𝒛 ‖2 𝑧𝑧 57

где 𝐴𝐴′∗ [𝜅𝜅(𝜉𝜉)] – сопряженный к A′ [𝜅𝜅(𝜉𝜉)] оператор, определяемый условием для -мерного вектора 𝒚𝒚 = {𝑦𝑦𝑞𝑞 , 𝑞𝑞 = 1 ÷ 𝑄𝑄}: 〈𝒚𝒚|A′ [𝜅𝜅(𝜉𝜉)]∆𝜅𝜅(𝜉𝜉)〉𝑅𝑅𝑄𝑄 = 〈𝐴𝐴′∗ [𝜅𝜅(𝜉𝜉)]𝒚𝒚|∆𝜅𝜅(𝜉𝜉)〉𝑋𝑋 Здесь 𝑋𝑋 – функциональное пространство для распределения приращения к коэффициенту пьезопроводности ∆𝜅𝜅(𝜉𝜉). В дискретном представлении для нашего случая это пространство 𝑅𝑅 С (где С – размерность ∆𝜅𝜅(𝑖𝑖, 𝑗𝑗) и Q – размерность 𝑦𝑦), а условие (3.3.10)(3.3.10) примет вид: �′ [𝜅𝜅(𝜉𝜉)]∆𝜅𝜅(𝜉𝜉)〉𝑅𝑅𝑄𝑄 = 〈𝐴𝐴̅′∗ [𝜅𝜅(𝜉𝜉)]𝒚𝒚|∆𝜅𝜅(𝜉𝜉)〉𝑅𝑅С 〈𝒚𝒚|A Подставляя в уравнение (3.3.11) выражение (3.3.5), получим: 𝑄𝑄 (𝑖𝑖,𝑗𝑗)=(𝑖𝑖,𝑗𝑗)𝑞𝑞𝑞𝑞 𝑞𝑞=1 (𝑖𝑖,𝑗𝑗)=(𝑖𝑖,𝑗𝑗)𝑞𝑞𝑞𝑞 � �𝑦𝑦𝑞𝑞 Для � того 𝑙𝑙(𝑖𝑖, 𝑗𝑗)∆𝜅𝜅(𝑖𝑖, 𝑗𝑗) �= 3𝜅𝜅 2 (𝑖𝑖, 𝑗𝑗) чтобы (𝑖𝑖,𝑗𝑗)=(𝑖𝑖,𝑗𝑗)𝑞𝑞𝑞𝑞 𝑄𝑄 � � �𝑦𝑦𝑞𝑞 (𝑖𝑖,𝑗𝑗)=(𝑖𝑖,𝑗𝑗)𝑞𝑞𝑞𝑞 𝑞𝑞=1 приведенное 𝑙𝑙(𝑖𝑖, 𝑗𝑗)∆𝜅𝜅(𝑖𝑖, 𝑗𝑗) � 3𝜅𝜅 2 (𝑖𝑖, 𝑗𝑗) равенство имело смысл следует специальным образом организовывать суммирование в правой части, что в конечном итоге определяет вид оператора 𝐴𝐴′̅ ∗ [𝜅𝜅(𝜉𝜉)]. Значение оператора 𝐴𝐴̅′∗ [𝜅𝜅(𝜉𝜉)] на -мерном векторе 𝒚𝒚 есть значения вектора 𝑔𝑔(𝑖𝑖, 𝑗𝑗) для индексов (𝑖𝑖, 𝑗𝑗) = {(𝑖𝑖, 𝑗𝑗)𝑞𝑞𝑚𝑚 ÷ (𝑖𝑖, 𝑗𝑗)𝑞𝑞𝑞𝑞 , 𝑞𝑞 = 1 ÷ 𝑄𝑄}. Для индексов (𝑖𝑖, 𝑗𝑗) выделяются значения 𝑞𝑞(𝑖𝑖, 𝑗𝑗), определяющие траектории 𝐿𝐿𝑞𝑞 , в интервалы которых попадают эти индексы. Далее конструируем выражение для вычисления значения сопряженного оператора на векторе 𝒚𝒚: 𝑄𝑄 𝑔𝑔(𝑖𝑖, 𝑗𝑗) = 𝐴𝐴′̅ ∗ [𝜅𝜅(𝜉𝜉)]𝒚𝒚 = � �𝑦𝑦𝑞𝑞(𝑖𝑖,𝑗𝑗) 𝑞𝑞(𝑖𝑖,𝑗𝑗) 𝑙𝑙𝑞𝑞(𝑖𝑖,𝑗𝑗) (𝑖𝑖, 𝑗𝑗) � 3𝜅𝜅 2 (𝑖𝑖, 𝑗𝑗) Здесь 𝑙𝑙𝑞𝑞(𝑖𝑖,𝑗𝑗) (𝑖𝑖, 𝑗𝑗) – это длина участка траектории 𝐿𝐿𝑞𝑞 от начальной точки до точки (𝑖𝑖, 𝑗𝑗), в которой выполняется расчет 𝑔𝑔(𝑖𝑖, 𝑗𝑗).Уравнение (3.3.13) при 𝒚𝒚 = 𝝋𝝋𝒛𝒛 примет вид: 58

𝛿𝛿𝐴𝐴𝑞𝑞1 [𝜅𝜅(𝜉𝜉)] 𝛿𝛿𝐴𝐴𝑄𝑄 [𝜅𝜅(𝜉𝜉)] ⋯ 𝑧𝑧 𝜑𝜑𝑞𝑞1 𝛿𝛿∆𝜅𝜅1𝑧𝑧 (𝑖𝑖, 𝑗𝑗) 𝛿𝛿∆𝜅𝜅1𝑧𝑧 (𝑖𝑖, 𝑗𝑗) ⎞ ⎛ 𝐴𝐴′∗ [𝜅𝜅(𝜉𝜉)]𝝋𝝋𝒛𝒛 = ⎜ ⋮ ⋱ ⋮ ⎟ × � ⋮𝑧𝑧 � 𝛿𝛿𝐴𝐴𝑞𝑞1 [𝜅𝜅(𝜉𝜉)] 𝛿𝛿𝐴𝐴𝑄𝑄 [𝜅𝜅(𝜉𝜉)] 𝜑𝜑𝑄𝑄 ⋯ 𝛿𝛿∆𝜅𝜅𝐶𝐶𝑧𝑧 (𝑖𝑖, 𝑗𝑗) ⎠ ⎝ 𝛿𝛿∆𝜅𝜅𝐶𝐶𝑧𝑧 (𝑖𝑖, 𝑗𝑗) После того как итерационный процесс завершился согласно условию в (3.3.7) и найдено решение ∆𝜅𝜅(𝑖𝑖, 𝑗𝑗) уравнения (3.3.5) для (3.3.6), служащее новым приближением для распределения пьезопроводности в узлах, определяемых траекториями вдоль сети, на области 𝑆𝑆 рассчитываются новые траектории 𝐿𝐿1𝑞𝑞 , на основе которых весь процесс повторяется. 3.4. Апробация вычислительной схемы гидродинамической томографии В данном разделе представлены результаты работы алгоритма вычислительной схемы гидродинамической томографии для различных моделей месторождений. Эксперимент № 1. Расчет интервального времени движения особых точек движения сигнала между скважинами на примере однородной среды нефтегазового месторождения. Начальные условия: сетка месторождения из двух скважин (см. рисунок 3.4.1), где скважина № 1 – является источником сигнала (давления) и № 2 – приемник сигнала. Карта среды с пьезопроводности 59 постоянным коэффициентом

Сетка месторождения для эксперимента № 1 Результаты расчета интервального времени для эксперимента № 1 Результаты расчета эксперимента представлены выше (см. рисунок 3.4.2). Траектория и интервальное время рассчитаны с использованием разработанного алгоритма на программном модуле. Эксперимент № 2. Расчет интервального времени движения особых точек движения сигнала между парами скважинам на примере неоднородной среды нефтегазового месторождения. 60

Начальные условия: сетка месторождения из 6-и скважин где все скважины последовательно являются источниками сигнала (давления) остальные 5 скважин – приемники сигнала. Карта среды в межскважинном пространстве имеет зоны с повышенным фильтрационным сопротивлением. Получены ожидаемые результаты (пути и времена), что свидетельствует о корректной и эффективной работе алгоритма (см. риунок 3.4.3). А Б В Г 61

Д Е Результаты работы алгоритма: А – скважина-источник № 1; Б – скважина-источник № 2; В – скважина-источник № 3; Г – скважина-источник № 4; Д – скважина-источник № 5; Е – скважина-источник № 6 Эксперимент № 3. Расчет интервального времени движения особых точек движения распределения сигнала между скважинами фильтрационного и сопротивления пространственного (коэффициент пьезопроводности) на примере тестовой модели нефтегазового месторождения. Начальные условия: сетка месторождения из 5 скважин, где каждая скважина последовательно является источником давления. Для нулевого приближения распределения фильтрационного сопротивления сгенерирована однородная среда 𝜅𝜅 0 (𝜉𝜉) (см. рисунок 3.4.4-А). Заданы тестовые наблюдаемые интервальные времена 𝝉𝝉� между парами скважин, 𝑄𝑄 = 20. На первой итерации по нулевому приближению находим траектории движения возмущения 𝑳𝑳𝒒𝒒 (см. рисунок 3.4.4-Б) между парами скважин и рассчитываем соответствующие этим траекториям интервальные времена 𝝉𝝉𝒒𝒒 . 62

А Б А – Сетка скважин и нулевое приближение среды; Б – Траектории движения сигнала между скважинами 𝐿𝐿𝑞𝑞 на первой итерации После завершения итерационного процесса и нахождения ∆𝜅𝜅(𝑖𝑖, 𝑗𝑗), определяем новые значения (Рисунок 3.4.5 - А). На рисунке 3.4.5 – Б представлен график динамики относительной погрешности 𝝋𝝋 между тестовыми и рассчитанными интервальными временами внутри итерационного процесса на первой итерации. На рисунках 3.4.6 – А и 3.4.6 – Б представлены результаты оптимизации среды после двух и десяти итераций соответственно. Б А А – Результаты работы оптимизационного процесса (7) на первой итерации; Б – Динамика средней относительной погрешности 𝝋𝝋 на первой итерации 63

А Б А - Результаты оптимизационного процесса на второй итерации; Б - Результат оптимизации среды после 10 итераций На рисунке 3.4.7 представлен график динамики средней относительной погрешности по ∆𝝉𝝉 после нахождения траекторий движения сигнала между скважинами. Далее происходит оценка этой величины и принятие решения о продолжении процесса вычислений. Динамика относительной погрешности по ∆𝝉𝝉 64

Результаты расчета пространственного распределения фильтрационного распределения вдоль траекторий движения сигнала между парами скважин для тестовой модели показал, что выбранный алгоритм оптимизации может быть применен для поиска 𝜅𝜅 = 𝜅𝜅 0 + ∆𝜅𝜅 на основе выбранного начального приближения 𝜅𝜅 0 . 3.5. Заключение Метод гидродинамической томографии проницаемых пластов нефтегазовых месторождений реализуется через создание депрессии в скважине, определение времени распространения нового установившегося уровня давления во всех скважинах гидродинамически связанных с источником депрессии и обработка рассчитанной матрицы времен по всем парам скважин. Времена распространения давления определяются с помощью принципов динамического программирования через поиск оптимальных траекторий. После необходимых измерений и получения матрицы интервальных времен возможна постановка задачи о построение карты пространственного распределения выраженного фильтрационного через коэффициент сопротивления проницаемого пьезопроводности, пласта, характеризующий процессы распространения давления внутри исследуемой среды. Задача решается через решение обратного оператора уравнения времени распространения экстремальной точки на кривой восстановления давления, между всеми парами скважин, использующихся в томографическом эксперименте. Рассматриваемый метод гидродинамической томографии позволяет проводить прогноз фильтрационных характеристик проницаемых пластов, контролировать участки возникновения низко проницаемых зон внутри продуктивных пластов нефтегазовых месторождений и устранять их. Своевременное устранение низко проницаемых зон позволит поддерживать объемы добываемых углеводородов на необходимом уровне, что особенно важно для месторождений, находящихся на поздних этапах разработки. 65

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

ГЛАВА 4. МЕТОДЫ ОЦЕНКИ РЕЗУЛЬТАТОВ РАБОТЫ ГИДРОДИНАМИЧЕСКОЙ ТОМОГРАФИИ Особенность гидродинамической данных используемых томографии при проницаемых решении пластов задачи приводит к возникновению существенных погрешностей в их значениях. Кроме того, выполнение алгоритма томографии реализуется с помощью решения интегрального уравнения для поля интервальных времен что является некорректно поставленной задачей. Поэтому существование погрешностей в данных приводит к значимым изменениям в вариантах решения обратной задачи гидродинамической томографии. Реальный интерес в этих условиях для практического применения рассчитываемого параметра, будут иметь для которых только те выполняется распределения оценка их погрешностей. Однако задача оценки погрешностей в решениях некорректных задач является логически недоопределенной и требует развития новых, специализированных подходов. Логическая недоопределенность состоит в том, что в любом диапазоне неопределенности данных можно найти элемент, для которого значение обратного оператора неограниченно. Из-за чего строго говорить о погрешности в решении, наследуемой погрешностями в данных невозможно без изменения принципов рассмотрения операторов. 4.1. Построение атрибута нечеткой меры по результатам гидродинамической томографии В разделе описана схема расчета атрибута нечеткой меры на основе интервальных оценок [108] фильтрационного сопротивления проницаемого плата, выраженного через коэффициент пьезопроводности по результатам работы метода гидродинамической томографии. распределение коэффициента пьезопроводности к (ξ ) Пространственное выражает сопротивление среды из-за которого зависит интервальное время движения особой точки (точки перегиба на КВД) восстановления давлений, между исследуемыми 67

скважинами. Перемещение особой точки происходит по траектории L ( rn , rend ) , соединяющей пару скважин, находящихся в координатах точек на сети ( rn , rend ) . Эта траектория такова, что интервальное время перемещения особой точки между ( rn ) и ( rend ) рассчитывается на следующему уравнению и должно быть минимально: τ ( rn, rend ) = Требованием минимальности l (ξ ) d ξ 3к (ξ ) L ( rn , rend ) ∫ (4.1.1) интервального времени определяется L ( rn , rend ) . l (ξ ) - длина пройденного пути от траектория движения особой точки L r ,r начальной точки rn до промежуточной ξ по траектории ( n end ) . Метод гидродинамической томографии состоит в поиске (приближенного) пространственного распределения коэффициента пьезопроводности к (ξ ) по сети, как функции пространственной координаты ξ по результатам расчетов всех = τ доступных для i τ ( rni , rend {τ= ). i } ;τ i измерения Эта интервальных задача решается времен пар традиционным скважин методом последовательных приближений на основе линеаризации (4.1.1) в окрестности начального приближения распределения коэффициента пьезопроводности. По определенному начальному приближению распределению значений параметра к 0 (ξ ) к устанавливаемому пьезопроводности, проводим линеаризацию уравнения (4.1.1) в окрестности начального приближения: { } ∆ τ =τ − τ = ∆ τ 0 0 0 i  l (ξ ) ∆к (ξ ) d ξ   ≈ − ∫  2 i 3  к 0 (ξ )    L ( rni ,rend ) (4.1.2) к= (ξ ) к 0 (ξ ) + ∆к (ξ ) Это система коэффициентов ∆к (ξ ) линейных алгебраических уравнений относительно , может быть выражена в матричной виде. Матричное представление уравнения (4.1.2) имеет вид: A(К 0 )∆К =∆ 0 τ 68 (4.1.3)

Метод последовательных приближений состоит в нахождении ∆К из (4.1.3), проверки качества результата 1 К= К 0 + ∆К и повторения при необходимости расчетов уже с новым начальным приближением параметров среды равным К1 , полученном на прошлой итерации процесса. Процесс длиться до получения такого распределения параметра К , при котором невязка между интервальными временами принимает удовлетворительное значение. Оператор, реализующий описанную схему в его матричном представлении будем записывать в форме: A−1 (К) τ = К (4.1.4) Отображение, вытекающее из (4.1.4) и предоставляющее возмоность i τ ( rni , rend {τ= ) прихода i } ;τ i = τ провести расчеты поля интервальных времен максимальных значений восстановления депрессии для системы исследуемых пар скважин (источник-приемник), которые зависят от распределения коэффициента пьезопроводности определяем: A(К) = τ , (4.1.5) Оператор A (К) τ = К является приближенным обратным [12] к «прямому» −1 оператору (4.1.5) и устанавливает способ решения (4.1.5) для определения матрицы среды коэффициента К пьезопроводности по экспериментальным (или полученным из математической модели эксплуатации = τ месторождения) интервальным временам Оператор компонентов (4.1.5) показывает матрицы К и i τ ( rni , rend {τ= ). i } ;τ i взаимосвязь параметров параметров матрицы среды из экспериментальных интервальных времен τ . Информационный граф, связывающий каждый элемент из К с каждым элементом из τ , оснащен весами для ребер, определяющими «коэффициент передачи» данных от элемента К ij A= граф оператора. Для линейных операторов роли весовых значений выступают к {a }; A ( К ) ≈ ∑ a элементы 69 τ i . Это информационный ijn ijn jn матрицы aijn . К ij= τ i , в Пример

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

построения оператора A и не являются интервалами доверия в общем случае. Изменение модельных представлений ведет к изменению конструкции оператора A , и как следствие, изменению вычисляемых интервалов доверия. В случае нелинейного оператора, к числу которых относится (4.1.4) отделение интервалов доверия от исходных данных неправомерно. Для поиска интервалов доверия к определяемым параметрам – коэффициентам пьезопроводности в рамках построенной модели, также следует рассчитать параметры с переносимыми данными для верхней границы интервала доверия по временам, далее повторить действия с данными времен, которые соответствуют нижней границе интервала доверия и рассмотреть разность их значений. Для исходных доверительных вычислениями экспериментальных интервалов Δτ = {∆τ i } данных (не строится τ путать с матрица промежуточными ∆ 0 τ в (2)) и рассчитываются два новых массива входных данных: τ sup = τ + Δτ = τ inf = τ − Δτ = {τ i + ∆τ i }; {τ i − ∆τ i }. (4.1.6) После по алгоритму, реализующему оператор гидродинамической томографии, находятся соответствующие матрицы параметров K sup и K inf с данными τ sup и τ inf . Интервал доверия ∆𝚱𝚱 рассчитывается как разница между двумя этими решениями: ∆K= В роли интервала доверия K sup − K inf могут (4.1.7) выступать соответствующие размерности данных атрибуты нечеткой меры [50] для исходных данных. Итогом вычислений будут служит атрибуты нечеткой меры для определяемых коэффициентов пьезопроводности. Рассмотрим в качестве тестового примера для апробации метода оценки результатов работы гидродинамической томографии экспериментальное месторождение (Рис. 4.1.2), на сетке которого расположены 28 скважин. Все скважины по очереди были использованы в качестве источника депрессии на 71

проницаемый пласт, и были смоделированы интервальные времена отклика на эту депрессию 𝝉𝝉 в соседних скважинах-приемниках. Всего в эксперименте по томграфическому исследованию пласта приняли участие 756 пар скважин, с целью охвата максимально возможной территории месторождения: 𝝉𝝉 = {𝜏𝜏𝑖𝑖 , i=1÷756}. Также следует отметить, что в эксперименте принято считать, все скважины гидродинамически связанными между собой, то есть депрессия, создаваемая в любой скважине сети, приводит к изменению режима работы во всех оставшихся. В качестве начального распределения фильтрационного сопротивления, выраженного через коэффициент пьезопроводности, избрана однородная среда (𝜅𝜅(𝜉𝜉(𝑥𝑥, 𝑦𝑦))=0,3), для которой представлена цветовая легенда по рассчитываемым значениям фильтрационного споротивления. Значения близкие к 1 соответствуют зонам аномального фильтрационного сопротивления (низкое значение коэффициента пьезопроводности), а значение близкие к нулю – соответствуют высокопроницаемым зонам. Сетка скважин экспериментального месторождения с заданным начальным приближением среды К На первом шаге вычислительного схемы алгоритма гидродинамической томографии находятся траектории движения волны депрессии (особой точки на кривой восстановления давления) между парами скважин 𝑳𝑳𝟎𝟎 и вычисляются 72

интервальные времена 𝝉𝝉𝟎𝟎 для начального приближения параметров среды К (Рис. 4.1.3). Траектории движения особой точки восстановления давления между парами скважин на первой итерации алгоритма гидродинамической томографии По итогам работы алгоритма на первой итерации гидродинамической томографии получаем новое итоговое пространственное распределение коэффициента пьезопроводности 𝑲𝑲𝟏𝟏 (Рис. 4.1.4), которое рассматривается в качестве начального приближения среды на следующем этапе вычислений. 73

Пространственное распределение коэффициента пьезопроводности по итогам первой итерации 𝑲𝑲𝟏𝟏 Финальный результат гидродинамической томографии продуктивного пласта пласта нефтегазового месторождения показан на рисунке 4.1.5. Заключение о завершении вычислительного процесса гидродинамической томографии принимается по результатам расчета отклонения модельных (экспериментальных) интервальных времен 𝝉𝝉 и времен 𝜏𝜏 𝑛𝑛 , полученных на итерации n. Итоговый результат гидродинамической томографии 74

После решения задачи гидродинамической томографии для интервальных времен 𝝉𝝉, строим новые массивы времен уже содержащие в себе доверительные 𝚫𝚫𝝉𝝉. Доверительный интервал для времен 𝝉𝝉 в эксперименте задан с уровнем доверия 0,8 от начальных интервальных времени. Поэтому разброс данных составлял 𝟎𝟎. 𝟐𝟐 от значений интервального времени в каждой используемой для расчетов паре скважин. Верхняя и нижняя границы доверительного интервала с уровнем доверия 0,8 𝝉𝝉𝒔𝒔𝒔𝒔𝒔𝒔 и 𝝉𝝉𝒊𝒊𝒊𝒊𝒊𝒊 рассчитывались следующим образом: τ sup = τ + (1 − 0,8) τ; (4.1.8) τ inf = τ − (1 − 0,8) τ. После для каждого массива времен 𝝉𝝉𝒔𝒔𝒔𝒔𝒔𝒔 и 𝝉𝝉𝒊𝒊𝒊𝒊𝒊𝒊 была заново решена задача гидродинамической томографии и рассчитаны соответствующие новых временам пространственные распределения коэффициента пьезопроводности 𝚱𝚱 𝒔𝒔𝒔𝒔𝒔𝒔 (Рис. 4.1.6) и 𝚱𝚱 𝒊𝒊𝒊𝒊𝒊𝒊 (Рис. 4.1.7). Пространственное распределение коэффициента пьезопроводности 𝚱𝚱 𝒔𝒔𝒔𝒔𝒔𝒔 для времен 𝝉𝝉𝒔𝒔𝒔𝒔𝒔𝒔 75

Пространственное распределение коэффициента пьезопроводности 𝚱𝚱 𝒊𝒊𝒊𝒊𝒊𝒊 для времен 𝝉𝝉𝒊𝒊𝒊𝒊𝒊𝒊 Интервал доверия 𝚫𝚫𝚫𝚫 (Рис. 4.1.8) для распределения коэффициента пьезопроводности 𝚱𝚱 между двумя решениями Κ 𝑠𝑠𝑠𝑠𝑠𝑠 и Κ 𝑖𝑖𝑖𝑖𝑖𝑖 рассчитывался по приведенному выше правилу (4.1.7). На краях сетки количество траекторий движения особой точки мало, в следствии чего в процессе расчетов эти участки мало изменялись относительно нулевого приближения, либо не изменялись вовсе. Именно с этим связаны пониженные значения интервала доверия. Это обстоятельство следует иметь в виду увязывая интервал доверия 𝚫𝚫𝚫𝚫 с плотностью сети траекторий (Рис. 4.1.9) в соответствующих зонах модели. Количество траекторий в ячейках сети максимального значения. 76 нормировано относительно

Интервал доверия 𝚫𝚫𝚫𝚫 коэффициента пьезопроводности 𝚱𝚱 Частота траекторий особой точки в ячейках сетки 4.2. Построение нечеткой модели распределения фильтрационного сопротивления проницаемого пласта Нечеткая величина [79, 80, 85, 86, 107] 𝜅𝜅 полностью определяется своей функцией принадлежности 0 ≤ 𝜇𝜇(𝜅𝜅) ≤ 1, которая имеет смысл меры достоверности, возможности того, что определение этой величины приведет к 77

значению 𝜅𝜅. В нашем случае, функция принадлежности 𝜇𝜇 для нечеткой величины коэффициента пьезопроводности 𝜅𝜅, рассчитывается в ячейках 𝜉𝜉(𝑖𝑖, 𝑗𝑗) сетки 𝑁𝑁 × 𝑀𝑀. Могут быть введены различные базовые функции принадлежности [61, 62, 68, 83, 87], комбинация которых образует итоговую функцию принадлежности [81, 82] для вариантов модели. Проведенные исследования [88, 90, 92], и результаты моделирования позволили в качестве наилучшей формы представления неопределенной информации для рассматриваемых задач выбрать экспоненциальную – гауссову функцию [84, 91]. Ее использование предполагает следующее. Во-первых, функция принадлежности должна учитывать изменение достоверности значений параметра 𝜅𝜅 в точке (𝑖𝑖, 𝑗𝑗), в которой рассчитывается значение возможности по мере удаления от значения 𝜅𝜅 в других точках сетки (𝑖𝑖2, 𝑗𝑗2): 𝜇𝜇1 �𝜅𝜅(𝑖𝑖, 𝑗𝑗)� = 𝑒𝑒𝑒𝑒𝑒𝑒 �− 2 �𝜅𝜅(𝑖𝑖2,𝑗𝑗2)−𝜅𝜅(𝑖𝑖,𝑗𝑗)� 2𝜎𝜎 2 �. (4.2.1) Во-вторых, необходимо учитывать отдаленность точки расчета от (𝑖𝑖2, 𝑗𝑗2): 𝜇𝜇2 �𝜅𝜅(𝑖𝑖, 𝑗𝑗)� = 𝑒𝑒𝑒𝑒𝑒𝑒 �− (𝑖𝑖2−𝑖𝑖)2 +(𝑗𝑗2−𝑗𝑗)2 2𝜎𝜎 2 �. (𝑖𝑖, 𝑗𝑗) (4.2.2) Этим обеспечивается уменьшение значения достоверности функции правдоподобия по мере удаления пробегаемых точек сетки от рассматриваемой. В формулах (4.2.1) и (4.2.2) параметр 𝜎𝜎 назовем эффективным, он характеризует степень разброса значений 𝜅𝜅. Для получения корректного результата необходимо рассматривать значения алгебраического произведения функций (1) и (2) отдельно по всем точкам сетки с последующим суммированием. Таким образом, итоговая функция принадлежности 𝜇𝜇, для нечеткой величины пьезопроводности 𝜅𝜅, вычисляется следующим образом: 1 𝜇𝜇�𝜅𝜅(𝑖𝑖, 𝑗𝑗)� = ∑𝑁𝑁,𝑀𝑀 𝑖𝑖2,𝑗𝑗2=1 �𝜎𝜎√2𝜋𝜋 ∙ 𝑒𝑒𝑒𝑒𝑒𝑒 �− 𝜅𝜅(𝑖𝑖2,𝑗𝑗2)≠0 Распределение месторождения �𝜅𝜅(𝑖𝑖2,𝑗𝑗2)−𝜅𝜅(𝑖𝑖,𝑗𝑗)� 2𝜎𝜎2 2 коэффициента находится с � ∙ 𝑒𝑒𝑒𝑒𝑒𝑒 �− (𝑖𝑖2−𝑖𝑖)2 +(𝑗𝑗2−𝑗𝑗)2 2𝜎𝜎2 пьезопроводности помощью 78 технологии коэффициента �� (4.2.3) исследуемого гидродинамической

томографии проницаемой среды согласно [4]. А соответствующие ему значения достоверности выраженной в функции принадлежности по формуле (3). Рассмотрим работы функции принадлежности [98] на примере модели нефтегазового месторождения, которая представлена в виде проекции поверхности исследуемого проницаемого пласта и разбита сетью ячеек размерностью 21 × 21 (рисунок 4.2.1–А). В пределах сети расположены пять скважин, для каждой пары скважин определены интервальные времена движения «сигнала» депрессии, от скважины-источника до скважин-приемников. В качестве нулевого приближения выбрана однородная среда (значение коэффициента пьезопроводности во всех ячейках 𝜅𝜅 = 0,3). На рисунке 4.2.1-А также представлены лучи движения реперной точки кривой восстановления давления между парами скважин на первой итерации гидродинамической томографии, вдоль которых были рассчитаны новые значения коэффициента пьезопроводности проницаемого пласта нефтегазового месторождения. Результат работы гидродинамической томографии представлен на рисунке 4.2.1-Б. А Б А – Сетка скважин с нулевым приближением среды; Б – Поле значений коэффициента пьезопроводности по результатам работы томографии без интерполяции 79

Поле данных, представленное на рисунке 4.2.1-Б служит исходными данными для построения функции принадлежности (3). На рис. 4.2.2 приведены поля достоверности пьезопроводности пространственного проницаемого пласта, распределения рассчитанные коэффициента по функции принадлежности (3) с разными значениями эффективного параметра 𝜎𝜎. 𝜎𝜎 = 0,5 𝜎𝜎 = 1,5 А Б 𝜎𝜎 = 1 𝜎𝜎 = 2,5 В Г Поле достоверности распределения коэффициента пьезопроводности при различных значениях эффективного параметра 80

Наиболее предпочтительные результаты расчета были получены со значением эффективного параметра, равным единице (Рисунок 4.2.2-В). Они хорошо отображают области с высоким и низким уровнями достоверности плотного и разреженного пространственного распределения коэффициента пьезопроводности. При меньшем значении эффективного параметра наблюдается очень резкий переход между высокими и низкими значениями достоверности (Рисунок 4.2.2-А). При значении эффективного параметра равном или более двух наблюдается практически одинаковая область распределения достоверности без учета плотности данных (Рисунок 4.2.2-Г). В качестве метода интерполяции для зон проницаемого пласта, в которых в ходе работы гидродинамической томографии не были определены значения коэффициента пьезопроводности, рассчитывается среднее по значениям в соседних ячейках. Поле значений коэффициента пьезопроводности с применением интерполяции На рисунке 4.2.3 приведено поле распределения коэффициента пьезопроводности после применения интерполяции для зон, в которых не были определены значения коэффициента пьезопроводности в результате работы гидродинамической томографии. 81

Следующий тестовый пример имеет другую конфигурацию расположения скважин (рисунок 4.2.4). В пределах сети расположены семь скважин, для каждой пары скважин определены интервальные времена движения «сигнала» депрессии, от скважины-источника приближения выбрана до скважин-приемников. однородная В качестве среда (значение нулевого коэффициента пьезопроводности во всех ячейках 𝜅𝜅 = 0,3). Результат работы гидродинамической томографии представлен на рисунке 4.2.5. Сетка скважин с нулевым приближением среды 82

Поле значений коэффициента пьезопроводности по результатам работы томографии без интерполяции Поле данных, представленное на рисунке 4.2.5, служит исходными данными для построения функции принадлежности. На рисунке 4.2.6 приведены поля достоверности пространственного распределения коэффициента пьезопроводности проницаемого пласта, рассчитанные по функции принадлежности (3) с разными значениями эффективного параметра 𝜎𝜎. 83

𝜎𝜎 = 0,5 𝜎𝜎 = 1,5 А Б 𝜎𝜎 = 1 𝜎𝜎 = 2,5 В Г Поле достоверности распределения коэффициента пьезопроводности при различных значениях эффективного параметра 84

Наиболее предпочтительные результаты расчета были получены со значением эффективного параметра, равным единице (Рисунок 4.2.6-В). Они хорошо отображают области с высоким и низким уровнями достоверности плотного и разреженного пространственного распределения коэффициента пьезопроводности. При меньшем значении эффективного параметра наблюдается очень резкий переход между высокими и низкими значениями достоверности (Рисунок 4.2.6-А). При значении эффективного параметра равном или более двух наблюдается практически одинаковая область распределения достоверности без учета плотности данных (Рисунок 4.2.6-Г). 4.3. Заключение Построение пространственного распределения эффективного фильтрационного сопротивления основано на построении пространственного распределения коэффициента пьезопроводности методом гидродинамической томографии. Для оценки интервала доверия к результатам расчетов в методе гидродинамической томографии исходными данными служат интервальные оценки уровня доверия к интервальным временам движения возмущений между парами скважин, служащие входными данными в методе гидродинамической томографии. Это соответствует третьему защищаемому положению: − Методы оценки достоверности результатов прогноза пространственного распределения фильтрационного сопротивления с применением теории нечетких множеств и с применением интервалов доверия исходных данных На основании интервальных оценок уровня доверия формируются пакеты данных нижней и верхней границ интервалов входных данных. Пользуясь информационным графом оператора гидродинамической принципами динамики нечетких мер для томографии и операторных уравнений, определяются нижние и верхние границы распределений коэффициента 85

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

ГЛАВА 5. ПРОГРАММНЫЙ КОМПЛЕКС МОДЕЛИРОВАНИЯ ПАССИВНОЙ ГИДРОДИНАМИЧЕСКОЙ ТОМОГРАФИИ ПРОНИЦАЕМОГО ПЛАСТА 5.1. Функциональные возможности Программный комплекс «Пассивная гидродинамическая томография» позволяет моделировать все этапы метода реконструкции фильтрационного сопротивления последовательно, как и отдельно каждый из этапов в зависимости от состава входных данных. На программный комплекс получено два свидетельства о государственной регистрации программы для ЭВМ. Функциональные возможности программного комплекса: 1. Импорт и экспорт данных. Работа с внешними данными в программе осуществляется с помощью excel шаблонов, которые содержат всю необходимую информацию. В программе возможен экспорт и импорт следующих данных: − параметры сетки карты месторождения; − координаты скважин; − история эксплуатации месторождения (скорости по дебиту и нагнетанию по скважинам); − оптимизационные параметры среды; − сетка пространственного распределения фильтрационного сопротивления, где каждая ячейка содержит свое значение. Размерность сетки сопротивления должна соответствовать размерности карты месторождения; − интервальные времена распространения депрессии от скважины- источника до скважин приемников − экспорт уровня достоверности пространственного распределения фильтрационного сопротивления. 87

− интервальные сопротивления по оценки для распределения данным томографической фильтрационного обработки данных гидродинамического прослушивания скважин. Пример шаблона для импорта данных по карте месторождения и сетки скважин 1. Прогноз работы скважин месторождения. Рассчитывается объем добычи по скважинам за интервал времени, установленный при загрузке истории эксплуатации скважин. Прогноз осуществляется после получения вектора параметров для которого наблюдается наименьшая невязка между модельными дебитами и дебитами из истории эксплуатации. Листинг модуля программного комплекса представлен в приложении Б.1 и Б.2. 2. Имитация депрессии в нагнетательной скважине и регистрация отклика в окружающих скважинах, что позволяет синтезировать интервальные времена прихода сигнала между скважинами 88 для метода гидродинамической

томографии. Листинг модуля программного комплекса представлен в приложении Б.3. 3. Поиск кратчайшего пути и расчет интервального времени движения сигнала между скважинами в неоднородной среде. 4. Реконструкция пространственного распределения фильтрационного проницаемого пласта нефтяного месторождения, где в качестве исходных данных используются: нулевое приближение распределения коэффициента пьезопроводности и интервальные времена, полученные в пункте 3 или загруженные из внешнего файла. 5. Построение нечеткой модели фильтрационного сопротивления проницаемого пласта по результатам работы гидродинамической томографии, где каждому значению фильтрационного сопротивления соответствует значение уровня достоверности расчетов. Листинг модуля программного комплекса представлен в приложении Б.4. 6. Построение интервальных оценок [111] фильтрационного сопротивления по данным томографической обработки данных гидродинамического прослушивания скважин. 5.2. Блок-схемы алгоритмов пассивной гидродинамической томографии В разделе представлены блок-схемы основных вычислтельных схем реализующихся в процессе решения задачи гидродинамичекой томграфии 89

Блок-схема алгоритма вычислительной схемы гидродинамической томографии 90

Начало Выбор нулевое приближение вектора ∆𝜅𝜅 0 (𝑖𝑖, 𝑗𝑗) Расчет сопряженного оператора Расчет параметра релаксации 𝛼𝛼 𝑧𝑧 Расчет значений вектора ∆𝜅𝜅 𝑧𝑧 (𝑖𝑖, 𝑗𝑗) Расчет приращения интервальных времен ∆𝝉𝝉𝒛𝒛𝒒𝒒 с учетом значений ∆𝜅𝜅 𝑧𝑧 (𝑖𝑖, 𝑗𝑗) Проверка значений 𝑧𝑧−1 (𝑖𝑖, 𝑗𝑗)� < 𝜀𝜀 �∆𝜅𝜅 (𝑖𝑖, 𝑗𝑗) − ∆𝜅𝜅 𝑧𝑧 Конец Блок схема внутреннего оптимизационного итерационного процесса 5.3. Пользовательский интерфейс программного комплекса На рисунке ниже приведена экранная форма главного окна программы, на котором отображается исследуемая среда (по координатная плоскость и сетка скважин. 91 умолчанию однородная),

Главное окно приложения (однородная среда нефтегазового месторождения и сетка скважин) На рисунке 5.3.2 продемонстрирована визуализация функции расчета оптимальных траекторий движения сигнала между скважинами. Главное окно приложения (траектории движения реперной точки кривой восстановления давления между парами скважин) 92

На рисунке 5.3.3 продемонстрировано отображение неоднородной среды [95], где более темные участки соответствуют слабопроницаемым зонам, а более светлые – сильнопроницаемые зоны. Главное окно приложения (неоднородная среда нефтегазового месторождения и сетка скважин) В дополнительном окне «Параметры карты месорождения» можно настроить основные параметры карты месторождения: ширина и высота сетки карты в реальных единицах, количество ячеек по осям X и Y и масштаб между ячейкой сетки и реальными единицами. 93

Форма настройки параметров карты месторождения В дополнительном окне «Координаты скважин» можно добавить \ удалить \ редактировать данные по скважинам. В случае загрузки данных сразу по нескольким скважинам используется шаблон excel-файла, который содержит информацию по каждой загружаемой скважине: № скважины, координаты скважины. 94

Форма настройки координат скважин месторождения На рисунке 5.3.6. Представлена окно просмотра истории эксплуатации скважин как в целом по месторождению так и отдельно по каждой скважине на протяжении загруженного интервала времени. В верхней полуплоскости отображается работа добывающих скважин (скорость добычи за единицу времени), а нижней полуплоскости – работа нагнетательных скважин (скорость закачки жидкости в пласт за единицу времени). 95

Окно просмотра истории эксплуатации по объемам добычи и нагнетания скважин/месторождения На рисунке 5.3.7 приведено окно приложения с параметрами среды по которым реализуется оптимизационного оптимизация процесса можно построенной импортировать модели. в Результаты excel-файл для последующего использования с следующих сеансах работы. Окно просмотра параметров модели эксплуатации месторождения На рисунке 5.3.8 представлено окно с результатами построения прогноза по работе скважин. Слева в таблице приведены следующие данные: модельный и реальный дебит, абсолютная и относительная погрешности, отдельные значения каждой компоненты модели (𝑄𝑄1 , 𝑄𝑄2 , 𝑄𝑄3 ). Справа на вкладке «Скважина» отображается графическое представление результатов прогноза. Верхний график отображает динамику изменения по компонентам дебита скважины. На нижнем графике приведено сравнение между дебитом из истории эксплуатации месторождения и дебитом рассчитанном с помощью модели на 96

выбранном временном интервале. На вкладке «Месторождение» представлены те же графике по всему месторождению без деления на скважины. На вкладке «Погрешность» можно отследить как изменяется погрешность между дебитами в различные интервалы времени [96]. Окно построения прогноза добычи по скважинам на заданном временном интервале На рисунке нагнетательной 5.3.9 скважине. представлено Выбирается окно имитации скважина, для депрессии которой в будет осуществляться имитация, устанавливается объем для нагнетания и его продолжительность. На графике фиксируются моменты времени, когда уровень добычи на скважинах-приемниках достигает максимального значения. 97

Окно моделирования интервальных времен распространения депрессии между скважинами 5.4. Заключение Для моделирования метода пассивной гидродинамической томографии проницаемого пласта программный комплекс нефтегазового в среде месторождения Microsoft Visual был Studio разработан на языке программирования C#. Разработанный программный комплекс реализует математическую модель эксплуатации скважин нефтегазового месторождения и вычислительную схему реконструкции пространственного распределения фильтрационного сопротивления проницаемого пласта. Это защищаемое положение номер четыре: − комплекс программ для реализации моделирования интервальных времен изменения по данным из истории эксплуатации скважин с целью прогнозирования пространственного распределения фильтрационного сопротивления и анализа результатов. 98

ЗАКЛЮЧЕНИЕ В ходе диссертационной работы выполнен литературный обзор методов гидродинамических исследований скважин, в частности гидродинамического прослушивания. Определены преимущества и недостатки рассмотренных методов и способы их решения с применением математического моделирования, численных методов и комплексов программ. Основные итоги выполненных в диссертации исследований состоят в следующем: Разработана нефтегазового имитационная математическая месторождения, гидродинамической томографии модель вычислительные для эксплуатации схемы моделирования пассивной данных с целью включения их в активную форму томографии. На математической модели реализуется имитация депрессии в одной скважине и регистрация отклика в окружающих скважинах, что позволяет синтезировать интервальные времена прихода сигнала между скважинами для метода гидродинамической томографии. Для сходимости итерационного процесса решаемой задачи используется параметр релаксации. Разработана вычислительная схема прогнозирования пространственного распределения фильтрационного сопротивления проницаемого пласта на основе данных пассивной гидродинамической томографии. Реализован метод поиска кратчайшего пути и расчета интервального времени на основе принципа Беллмана. Разработаны два способа оценки результатов прогноза пространственного распределения фильтрационного сопротивления с применением теории нечетких множеств, пространственного построена распределения функция коэффициента принадлежности для пьезопроводности, где каждому значению рассчитывается уровень достоверности. Второй способ основан на расчете интервальных оценок для коэффициента пьезопроводности по данным томографической обработки для данных гидродинамичекого прослушивания скважин исследуемого нефтегазового месторождения. 99

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

СПИСОК ЛИТЕРАТУРЫ 1. Вольпин С. Г. Анализ применения технологий гидродинамического исследования скважин в информационном обеспечении проектирования разработки [Текст] / Вольпин С. Г., Мясников Ю. А., Свалов А. В. // Нефтяное хозяйство. – 2002. – № 10. – С. 61–65. 2. Вольпин С. Г. Состояние гидродинамических исследований скважин в нефтедобывающей отрасли России [Текст] / Вольпин С. Г. // Нефтяное хозяйство. – 2003. – № 6. – С. 66–68. 3. Чодри А. Гидродинамические исследования скважин [Текст]: пер. с англ. В. А. Юдина, О. В. Ломакиной / Чодри А. – Москва.: Премиум инжиниринг, 2011. – 687 с. 4. Патент 2092691 РФ. Способ контроля фильтрационных потоков, формирующихся при разработке нефтяных месторождений со слоистонеоднородными пластами [Текст] / Кондаратцев С. А., Хасанов М. М., Хатмуллин И. Ф., [и другие]: Внедренческий научно-исследовательский инженерный центр. – № 95101668/03; заявление: 10.02.1995; публикация: 10.10.1997. 5. Патент 2229020 РФ. Способ выявления непроводящих элементов нефтяной залежи при ее эксплуатации [Текст] / Щацкий А. В., Колесов В. В., Щацкий Д. А. [и другие]: ЗАО Пангея. - № 2002129342; заявление: 05.11.2002; публикация: 20.05.2004. 6. Патент 2298647 РФ. Способ исследования нефтяных пластов [Текст] / Щацкий А. В., Колесов В. В., Щацкий Д. А. [и другие]: ЗАО Пангея. – № 2005111998; заявление. 22.04.2005; публикация: 10.05.2007. 7. Распопов Р. В. Регуляризация оценок гидродинамических параметров нефтеносного коллектора в технологиях группового гидропрослушивания: автореферат диссертации : 05.13.18 / Распопов Роман Владимирович. – Тюмень: ТГНУ, 2015. – 19 с. 101

8. Басниев К. С. Подземная гидромеханика [текст] / Басниев К. С., Дмитриев Н. М., Каневская Р. Д. [и другие]. – М.: Институт КИ. – 2006. – 487 с. 9. Каневская Р. Д. Математическое моделирование гидродинамических процессов разработки месторождений углеводородов [текст] / Каневская Р. Д. – М.: Институт КИ. – 2003. – 128 с. 10. Богачев К. Ю. Эффективное решение задач фильтрации вязкой сжимаемой многофазной многокомпонентной смеси на параллельных ЭВМ: дис. док. физ.-мат. наук: 05.13.18 / Богачев Кирилл Владимирович. – М., 2012. – 201 с. 11. Юдин Е. В. Применение математической модели для уточнения некоторых геологических и гидродинамических характеристик нефтяной залежи с использованием промысловых данных [текст] / Юдин. Е. В. – М: МФТИ., 2006. 12. Бакиров Э. А. Геология нефти и газа [учебник для вузов] / Под ред. Э. А. Бакирова, Ермолкин В. И., Ларин В. И. – М.: Недра, 1990. – 240 c. 13. Звездин В. Г. Нефтепромысловая геология [Учебно-методическое пособие] / Звездин В. Г. Перм: ПГУ, 2007. — 115 c. 14. Гавич И. К. Основы гидрогеологии. Гидрогеодинамика / И. К. Гавич, И. С. Зекцер, В. С Ковалевский [и др.]. – Новосибирск: Наука, 1983. – 242 с. 15. Мищенко И. Т. Скважинная добыча нефти [текст] // Мищенко И. Т.– Москва: Нефть и газ, 2003. – 806 с. 16. Баренблатт Г. И., Желтов Ю. П. Об основных уравнениях фильтрации однородных жидкостей в трещиноватых породах [текст] // ДАН СССР. – Т. 132. – №3. –1960. – С. 545–548 17. Борисов Ю. П. Определение параметров пласта при исследовании скважин на неустановившихся режимах с учетом продолжающегося притока жидкости [текст] // ВНИИ. – Выпуск 19. – 1959. –C. 115–133. 18. Борисов Ю. П., Баренблатт Г. И. Об определении параметров пласта по данным о восстановлении давления в остановленных скважинах [текст] // Борисов Ю. П., Баренблатт Г. И. [и другие]. –Известия АН СССР. – 1957. – № 11. –C. 84–91. 102

19. Борисов Ю. П., Блох С. С., Митюшо В. Н. Анализ некоторых методов обработки кривых восстановления давления в неоднородных пластах [текст] // ВНИИ. – 1970. – Выпуск 55. – C. 174–188. 20. Борисов Ю. П., Орлов B. C. Интерпретация данных восстановления давления и их использование при построении карт изобар [текст] // журнал «Нефтяное хозяйство». – 1957. – № 7. – C. 39–43. 21. Борисов Ю. П., Орлов Ф. Ф. Исследование нагнетательных скважин месторождения Контур-Так [текст] / Борисов Ю. П., Орлов Ф. Ф. // Нефтепромысловое дело. – 1979. – № 10. – C. 17–20. 22. Борисов Ю. П., Требин Ф. А., Мухарский Э. Д. определению параметров пласта по кривым восстановления давления [текст] // журнал «Нефтяное хозяйство». – 1958. – № 8. – C.38–45. 23. Щелкачев В. Н. Основы и приложения теории неустановившейся фильтрации. Часть 1 [монография] // Щелкачев В. Н. – Москва.: Нефть и газ, 1995. – 586 с. 24. Щелкачев В. Н. Разработка нефтегазоносных пластов при упругом режиме [текст] // Щелкачев В. Н. – М.: Гостоптехиздат, 1959. – 467 с. 25. Пикалов В. В., Преображенский Н. Г. Вычислительная томография и физический эксперимент // Успехи физических наук. – 1983. –Том 141-3. – С. 469-493 26. Кобрунов А. И. Теоретические основы гидродинамической томографии // Геофизический журнал. – 2015. – Выпуск 2. – С. 27–34. 27. Кобрунов А. И. Математическая модель томографии на давлениях при контроле за разработкой нефтяных месторождений // Известия Коми научного центра Уро РАН. – 2012. – Выпуск 4-12. – С. 82–86. 28. Марчук Г.И. Методы вычислительной математики. – Москва: Наука, 1980. – 534 с. 29. Болтянский В. Г. Математические методы оптимального управления. – Москва: Наука, 1969. – 408 с. 103

30. Енцов И. И. Определение гидродинамических параметров по данным исследования скважин на приток при неустановившемся режиме [Учебнометодическое пособие] // Енцов И. И. – Ухта: Ухтинский государственный технический университтет, 2010. – 11 с. 31. Лаврентьев М. М. Некорректные задачи математической физики и анализа / М. М. Лаврентьев, В. Г. Романов, С. П. Шишатский. – М.: Наука, 1980. – 286 с. 32. Воронин С. Г., Курносов Д. А., Корабельников М. И. Математическое моделирование эксплуатационной скважины в процессе оптимизации нефтедобычи [текст] Воронин С. Г., Курносов Д. А., Корабельников М. И. [и другие] // Вестник ЮУрГУ: серия Энергетика. 2005. №9 (49). URL:https://cyberleninka.ru/article/n/matematicheskoe-modelirovanieekspluatatsionnoy-skvazhiny-v-protsesse-optimizatsii-neftedobychi 33. Краснов В. А., Иванов В. А., Хасанов М. М. Помехоустойчивый метод оценки связности пласта по данным эксплуатации месторождений // Российская техническая нефтегазовая конференция и выставка SPE по разведке и добыче (Москва, 16–18 октября 2012 года). – SPE. – 162053. 34. Jong S. Kim, Larry W. Lake, Thomas F. Edgar. Integrated CapacitanceResistance Model for Characterizing Waterflooded Reservoirs // Proceedings of the 2012 IFAC Workshop on Automatic Control in Offshore Oil and Gas Production (May 31-June 1, 2012).- Norwegian University of Science and Technology, Trondheim.- 2012.- P. 19-24. 35 Valko P. P., Doublet L. E., Blasingame T. A. Development and application of the Multiwell Productivity Index (MPI), SPE, Texas A\& M U, March 2000. 36. De Sant Anna Pizarro, J. O. 1998. Estimating Injectivity and Lateral Autocorrelation in Heterogeneous Media. PhD dissertation. U. Of Texas, Austin, Texas. 37. Пархоменко Я. В. Мониторинг продуктивности многоскважинной системы при разработке месторождений углеводородов // Пархоменко Я. В. – М.: МФТИ. 104

38. Ozkan E. Performance of Horizontal Wells, PhD dissertation, U. of Tusla, Tusla, Oklahoma, 1988. 39. Бузинов С. Н., Умрихин И. Д. Влияние неоднородности пласта по напластованию на определение его параметров по данным наблюдения его нестационарной фильтрации [текст] // Ежегодник ВНИИ. – Москва, 1966. – C.307-321. 40. Бузинов С. Н., Умрихин И. Д. Исследование нефтяных и газовых скважин и пластов [текст] // Бузинов С. Н., Умрихин И. Д. – Москва: Недра, 1984. – 269 с. 41. Бузинов С. Н., Умрихин И. Д. Исследование пластов и скважин при упругом режиме фильтрации [текст] // Бузинов С. Н., Умрихин И. Д.– Москва: Недра, 1964. – 273 с. 42. Бузинов С. Н., Умрихин И. Д. Гидродинамические методы исследования скважин и пластов [текст] // Бузинов С. Н., Умрихин И. Д.. – Москва: Недра, 1973. – 246 с. 43. Умрихин И. Д., Днепровская Н. И. Состояние и проблемы гидродинамических исследований [текст] // Умрихин И.Д., Днепровская Н.И. [и другие] // Нефтяное хозяйство. – 1993. –№3. –С.55-57. 44. Умрихин И. Д., Днепровская Н. И. Исследование нефтяных скважин на нескольких режимах [текст] // Умрихин И. Д., Днепровская Н. И. [и другие]. – журнал «Нефтяное хозяйство». –1988. – № 7. – С.37-39. 45. Басниев К. С. Нефтегазовая гидромеханика [Учебное пособие для вузов] // Басниев К. С., Дмитриев, Н. М., Розенберг, Г. Д. – Москва: Институт компьютерных исследований, 205. - 544 с. 46. Кунцев В. Е., Кобрунов А. И., Мотрюк Е. Н. Технология оценки связности скважин на основе модели эксплуатации месторождения [текст] Кунцев В. Е., Кобрунов А. И., Мотрюк Е. Н. // Фундаментальные исследования. – 2015. – № 6–3. – С. 452–456. 47. Басниев К. С. Подземная гидравлика [Учебник для вузов] // Басниев К. С., Власов, А. М., Кочина [и другие]. -Москва: изд-во: Недра, 1986. - 303 с. 105

48. Кунцев В. Е., Кобрунов А. И., Мотрюк Е. Н. Вычислительная схема гидродинамической томографии [текст] Кунцев В. Е., Кобрунов А. И., Мотрюк Е. Н. // Фундаментальные исследования. – 2016. – № 7-2. – С. 230-235. 49. Тихонов, A. Н., Арсенин, В. Я. Методы решения некорректных задач [Учебник для вузов] // Тихонов, A. Н., Арсенин, В. Я. - Москва.: изд-во Наука, 1979. – 285 c. 50. Кобрунов А. И. Нечетко-атрибутный анализ решений операторных уравнений [текст] // А. И. Кобрунов. - Известия Самарского научного центра РАН. – 2017. – том 19, № 1(2). – с. 410–413. 51. Соловьев И. Г., Распопов Р.В. Регуляризация оценок параметров нефтяных коллекторов по условиям симметрии [текст] // Автоматизация, телемеханизация и связь в нефтяной промышленности. – 2010. – № 11. – С. 28-33. 52. Соловьев И. Г., Распопов Р.В. Устойчивое оценивание параметров коллекторов на основе ортогонализации [текст] // Вестник Кибернетики. – 2010. – № 9. – С. 20-27. 53. Стрекалов А. В. Математические модели гидравлических систем для управления системами поддержания пластового давления [текст] // Стрекалов А. В. - Тюмень: Дом печати, 2007. – 661 с. 129. 54. Тихонов A. Н., Гончарский А. В. Численные методы решения некорректных задач [Учебник для вузов] // Тихонов A. Н., Гончарский А. В., Степанов В. В., [и другие].- Москва.: Книга по требованию, 2012. – 228 c. 55. Форсайт Д., Молер К. Численное решение систем линейных алгебраических уравнений [текст] Форсайт Д., Молер К. – Москва.: Мир, 1969. – 168 с. 56. Хасанов М. М, Краснов В. А. Решение задачи о взаимодействии пласта со скважиной в условиях нестационарного притока [текст] Хасанов М. М, Краснов В. А., [и другие] // Научно-технический вестник «Роснефть».- 2007.- № 2.- С. 41-46. 57. Хатмуллин И. Ф., Латыпов А. Р. Принципы построения адаптивной постоянно действующей модели нефтяной залежи [текст] Хатмуллин И. Ф., 106

Латыпов А. Р. [и другие].- журнал «Нефтяное хозяйство». – 2004. – № 10. – С. 58-61. 58. Чарный И. А. Подземная гидрогазодинамика [текст] // Чарный И.А.- Москва: Современные нефтегазовые технологии, 2007.-436 с. 59. Щелкачев В. Н., Лапук Б. Б. Подземная гидравлика [учебник для вузов] // Щелкачев В. Н., Лапук Б. Б.– Ижевск: НИЦ Регулярная и хаотическая динамика, 2001.- 736 с. 60. Щуров В. И. Технология и техника добычи нефти [учебник для вузов] // Щуров В. И.- Москва: Недра, 1983.- 510 с. 61. Семухин М. В. Нечеткие оценки запасов нефти [текст] // Сборник докладов международной конференции по мягким вычислениям.- Санкт-Петербург, 2003.- Том 2.- С. 164-167. 62. Ягер Р. Нечеткие множества и теория возможностей. Последние достижения [текст] // перевод Кузьмин В. Е., Травник С. Н.- Москва: Радио и связь, 1986.- 408 с. 63. Сургучев М. Л. О нестационарных режимах заводнения нефтяных пластов [текст] // Сургучев М. Л., Цынкова О. Э.- журнал «Нефтяное хозяйство».- 1983.–-№7.- С. 26-29. 64. Терещенко С. А. Методы вычислительной томографии [текст] // Терещенко С.А.- М.: ФизМатЛит, 2004.- 319 с. 65. Маргот Г. Джеристон и Льюис Дж. Дурловски. Моделирование жидкости в нефтяных пластах. 66. Япа П., Женг Л. И Наката К. Моделирование подводных нефтегазовых струй и шлейфов. 67. Шечтер Д. С., Гуо Б. Математическое моделирование тяжести дренаж после закачки газа в трещиноватые коллекторы. 68. Zadeh L. Fuzzy sets / Information and Control. 1965, volume 8. Page 338-353. 69. Sarma P. Efficient real-time reservoir management using adjoint on base optimal control and model updating / Computational Geosciences Journal. 2006. – Volume 10. page 3-36. 107

70. Doren J. Reduced-order optimal control of water flooding using proper orthogonal decomposition / Computational Geosciences Journal. 2006. Volume. 10.- page 137-158. 71. Sturm W. Dynamic Reservoir Well Interaction / SPE Annual Technical Conference and Exhibition. Texas, 2004. 72. Echeverria D. A. Robust Scheme for Spatio-Temporal Inverse Modeling of Oil Reservoirs / 18th World IMACS / MODSIM Congress, Cairns, Australia 13-17 July 2009. page. 4206–4212. 73. Gao C. Literature Review on Smart Well Technology / Production and Operations Symposium, 31 March-3 April 2007, Oklahoma City, Oklahoma, USA. SPE, 2007. 74. Gao G. A Stochastic Optimization Algorithm for Automatic History Matching / G. Gao, G. Li, A. Reynolds // SPE Journal. 2007. Volume 12. Page 196-208. 75. Glandt C. Reservoir Aspects of Smart Wells / C. Glandt // SPE Latin American and Caribbean Petroleum Engineering Conference, 27-30 April 2003, Port-of Spain, Trinidad and Tobago. SPE, 2003. 76. Oliver D. Recent progress on reservoir history matching / Computational Geosciences. 2010. Volume. 15. № 1. P. 185-221. 77. Pashali A. Real Time Optimisation Approach for wells / Intelligent Energy Conference and Exhibition, 25-27 February, Amsterdam, The Netherlands. SPE, 2008. page 93-98. 78. Mamdani E. Application of fuzzy algorithms for control of simple dynamic plant / Proceedings of the IEE, 1974. Volume 121. Page 1585-1588. 79. Алексеев А. В. Применение нечеткой математики в задачах принятия решений [текст] / Алексеев А.В. // Методы и системы принятия решений.- Рига: РПИ, 1983.- C. 38-42. 80. Дорогобед А. Н. Геомоделирование в условиях неопределенности для задач газопромысловой отрасли: диссертация. канд. техн. Дорогобед Алена Николаевна. – Петрозаводск, 2016. – 152 с. 108 наук: 05.13.18 /

81. Алтунин Е. А. Нечеткие методы идентификации и управления процессами нефтегазодобычии: диссертация канд. техн. наук.: 05.13.01 // Алтунин Евгений Александрович.- Тюмень, 2002.- 203 с. 82. Алтунин А. Е.Методы определения функций принадлежности в теории размытых множеств [текст] Алтунин А. Е. // Труды ЗапсибНИГНИ.- Тюмень, 1980.- Выпуск154.- С. 62-72. 83. Алтунин А. Е., Семухин М. В. Применение теории нечеткости для оценивания технологических параметров [текст] // Алтунин А. Е., Семухин М. В. [и другие].- Проблемы нефти и газа. Тюмень, 1983. Выпуск 58.- С. 57-59. 84. Алтунин А. Е., Семухин М. В. Методические рекомендации по применению теории нечеткости в процессах контроля и управления объектами газоснабжения [текст] // Алтунин А. Е., Семухин М. В. [и другие].- Тюмень, 1983.- 136 с. 85. Заде Л. А. Новое в зарубежной науке. Понятие лингвистической переменной и его применение к принятию приближенных решений [текст] // Заде Л. А.- Москва: Мир, 1976.- 165 с. 86. Заде Л. А. Размытые множества и их применение в распознавании образов и кластер-анализе [текст] // Заде Л. А.- Классификация и кластер.- Москва, 1980 .- С. 208-247. 87. Кандель А. Нечеткие множества, нечеткая алгебра, нечеткая статистика [текст] // Кандель А., Байатт У. Труды американского общества инженеров радиоэлектроников.- 1978.-Том 66-12.- С.37-61. 88. Кобрунов А. И. Методы нечеткого моделирования при изучении взаимосвязей между геофизическими параметрами [текст] Кобрунов А. И. // Геофизика.- 2010.-Том 2.- C.17-23. 89. Кобрунов А. И. Метод изучения пространственного распределения фильтрационного сопротивления при эксплуатации нефтяных месторождений [текст] / Кобрунов А. И., Художилова А. Н. [и другие] // Нефтяное хозяйство.- 2013.- Том 8.- С. 58-60. 109

90. Кобрунов А. И. Метод нечетких петрофизических композиций при прогнозировании петрофизических параметров / Кобрунов А. И., Художилова А.Н., [и другие] // Вестник института геологии Коми научного центра УРО РОАН.- 2011.- №9.- С. 18-24. 91. Норвич A. M. Фундаментальное измерение нечеткости [текст] // Норвич A.M.– Москва: Радио и связь, 1986. – С.54-64. 92. Кобрунов А. И. Адаптация метода нечѐтких петрофизических композиций для определения подсчѐтных параметров Низевого месторождения // Кобрунов А. И., [и другие] // Электронный научный журнал «Нефтегазовое дело».- 2011.- №6.- С. 307-315. 93. Семухин М. В. Нечеткие оценки запасов нефти / Семухин М. В.// Сборник докладов международной конференции по мягким вычислениям. – С-Пб, 2003. – Том 2. – С. 164-167. 94. Семухин М. В. Теория нечетких множеств [Учебно-методическое пособие] // Семухин М. В.- Тюмень: ТюмГУ, 1999. – 50 с. 95. Кунцев В. Е. Пассивная гидродинамическая томография проницаемого пласта [Текст] / В. Е. Кунцев, А. И. Кобрунов, Е. Н. Мотрюк. – Свидетельство о государственной регистрации программы для ЭВМ № 2017662707; заявл. 24.07.2017; опубл. 15.11.2017. -1 с. 96. Кунцев В. Е. Прогнозирование работы скважин нефтегазового месторождения [Текст] / В. Е. Кунцев, А. И. Кобрунов, Е. Н. Мотрюк. – Свидетельство о государственной регистрации программы для ЭВМ № 2017619300; заявл. 12.04.2017; опубл. 21.08.2017. -1 с. 97. Кобрунов А.И., Кунцев В.Е., Мотрюк Е.Н. Методика и технология решения задачи пассивной гидродинамической томографии// Вопросы теории и практики геологической интерпретации геофизических полей: материалы 43-й сессии Международного семинара им. Д.Г. Успенского (г. Воронеж, 26-30 января 2016г.), Воронеж: Издательско-полиграфический центр «Научная книга», 2016.– С. 107-110 110

98. Кобрунов А. И., Мотрюк Е. Н., Кунцев В. Е. Построение нечеткой геологической модели на основе результатов гидродинамической томографии // Успехи современной науки. – 2016. – № 12–5. – С. 43–46 99. Кобрунов А. И., Кунцев В. Е., Мотрюк Е. Н. Технология оценки связности скважин на основе модели эксплуатации месторождения // Фундаментальные исследования. – 2015. – № 6–3. – С. 452–456 100. Кобрунов А. И., Кунцев В. Е., Мотрюк Е. Н. Вычислительная схема гидродинамической томографии // Фундаментальные исследования. – 2016. – № 7–2. – С. 230–235 101. Кобрунов А. И., Мотрюк Е. Н., Кунцев В. Е. Алгоритм поиска кратчайшего пути и интервального времени между скважинами на основе принципа Беллмана // Современные наукоемкие технологии. – 2016. – № 8–1. – С. 51–55 102. Кунцев В. Е. Математическая модель оценки связности скважин для метода гидродинамической томографии // Международный студенческий научный вестник. -2015. - №-3-3. – С 397–399. 103. Кобрунов А. И., Кунцев В. Е., Мотрюк Е. Н. Теоретические предпосылки и принципы реализации пассивной гидродинамической томографии проницаемого пласта // Ресурсы Европейского Севера. Технологии и экономика освоения. – 2015. – № 1. – С. 19–26. 104. Кобрунов А. И., Куделин С. Г., Художилова А. Н., Кунцев. Методы и технологии гидродинамической томографии проницаемых пластов // Вопросы теории и практики геологической интерпретации гравитационных, магнитных и электрический полей: материалы 42-й сессии Международного семинара им. Д. Г. Успенского. – Горный ин-т УрО РАН. – Пермь, 2015. – С. 105–107 105. Кунцев В. Е., Мотрюк Е. Н., Кобрунов А. И. Математическая модель пассивной гидродинамической томографии проницаемого пласта // Рассохинские чтения [Текст]: материалы международного семинара (4–5 февраля 2016 г). – В 2 ч. Ч. 2. – Ухта: УГТУ, 2016. – С. 128–133 106. Кунцев В. Е., Мотрюк Е. Н., Кобрунов А. И. Математическая модель оценки связности скважин по данным истории эксплуатации месторождения // 111

Рассохинские чтения [Текст]: материалы международного семинара (5–6 февраля 2015 г). – В 2 ч. Ч. 2. – Ухта: УГТУ, 2015. – С. 143–149 107. Кунцев В. Е., Мотрюк Е. Н., Кобрунов А. И. Построение нечетких моделей фильтрационного сопротивления проницаемого пласта // Рассохинские чтения [Текст]: материалы международной конференции (2-3 февраля 2017 г.). В 2 ч. Ч. 2 / под ред. Н. Д. Цхадая. - Ухта: УГТУ, 2017. - С. 117-120 108. Кобрунов А. И., Кунцев В. Е. Интервальные оценки для коэффициента пьезопроводности по данным томографической обработки данных гидродинамического прослушивания // Известия Самарского научного центра РАН. – 2017. - том 19, № 4. – С. 153-160. 109. Кунцев В. Е. Метод пассивной гидродинамической томографии проницаемого пласта // Физико–математическое моделирование систем [Текст]: материалы XVIII Международного семинара (г. Воронеж, 30 июня 2017 г.). Ч. 1. – Воронеж: ВГТУ, 2017. – С. 139–146. 110. Патент 2501948 РФ. Способ гидродинамической томографии проницаемости пласта // Кобрунов А. И. заявитель и патентообладатель ФГБО ВПО «УГТУ».– № 2012120976/03; заявление. 22.05.2012; опубл. 20.12.2013. 111. Кунцев В. Е. Построение интервальных оценок гидродинамической томографии проницаемого пласта [Текст] / В. Е. Кунцев, А. И. Кобрунов, Е. Н. Мотрюк. – Свидетельство о государственной регистрации программы для ЭВМ № 2018617703; заявление. 30.05.2018; опубл. 28.06.2018. -1 с. 112

ПРИЛОЖЕНИЯ 113

ПРИЛОЖЕНИЕ А Свидетельсва о регистрации и акты о внедрении Свидетельство о государственной регистрации программы для ЭВМ «Прогнозирование работы скважин нефтегазового месторождения» 114

Свидетельство о государственной регистрации программы для ЭВМ «Пассивная гидродинамическая томография проницаемого пласта» 115

Свидетельство о государственной регистрации программы для ЭВМ «Построение интервальных оценок гидродинамической проницаемого пласта» 116 томографии

Акты о внедрении 117

118

119

120

ПРИЛОЖЕНИЕ Б Текст программы «Пассивная гидродинамическая томография проницаемого пласта» Приложение Б.1. Программная реализация создания объекта нефтегазового месторождения с целью последующих расчетов параметров среды. class DevField { public MapField mapField; // карта нефтегазового месторождения public int countWells; // количество скважин месторождения public int countInjWells; public int countProdWells; public List<Well> listWells; // список скважин public List<PairWells> listPairWells; // Конструктор класса нефтегазового месторождения public DevField() { mapField = new MapField(); listWells = new List<Well>(); } // Функция добавления скважины в список public void AddWell(int prmNumberWell, int prmCoordX, int prmCoordY) { listWells.Add(new Well(prmNumberWell, prmCoordX, prmCoordY)); countWells++; } // Функция добвления нескольких скважин public void AddWells(List<Well> prmListWells) { for (int i = 0; i < prmListWells.Count; i++) { listWells.Add(prmListWells[i]); } } // Функция возвращает скважину по ее номеру public Well GetWell(int prmNumberWell) { return listWells.Find(x => x.numberWell == prmNumberWell); } // Процедура рандомной генерации двух скважин public void InitRandomPairWell() { listWells.Clear(); Random rand = new Random(); AddWell(1, rand.Next(1, mapField.sizeGridX), rand.Next(1, mapField.sizeGridY)); AddWell(2, rand.Next(1, mapField.sizeGridX), rand.Next(1, mapField.sizeGridY)); countWells = 2; } // Процедура генерации двух скважин по координатам public void InitUserPairWell(int prmX1, int prmY1, int prmX2, int prmY2, bool prmFlag) { if (prmFlag) listWells.Clear(); // чистим список скважин AddWell(1, prmX1, prmY1); 121

AddWell(2, prmX2, prmY2); countWells = 2; } // Функция проверки наличия скважины ко координатам ячейки public bool CheckCellWell(int prmX, int prmY) { if (listWells.Find(s => s.coordXWell == prmX && s.coordYWell == prmY) != null) return true; else return false; } // Функция возвращает номер скважины по координатам public int GetNumberWell(int prmCoordX, int prmCoordY) { return listWells.Find(x => x.coordXWell == prmCoordX && x.coordYWell == prmCoordY).numberWell; } public void DeleteWellAll() { listWells.Clear(); countWells = 0; } // Работа с параметрами public Model.ListParameters listDistanceWells; // список параметров: расстояние между скважинами // Процедура иниуиализации списков хранения параметров модели public void InitList() { listDistanceWells = new Model.ListParameters(); } // public void CalcListDistanceWells() { for (int w1 = 0; w1 < countWells; w1++) { Well well1 = listWells[w1]; for (int w2 = 0; w2 < countWells; w2++) { Well well2 = listWells[w2]; if (well1.numberWell != well2.numberWell) { listDistanceWells.AddParameter(well1.numberWell, well2.numberWell, CalcDistance(well1.coordXWell, well1.coordYWell, well2.coordXWell, well2.coordYWell)); } } } } // private double CalcDistance(int prmCoordX1, int prmCoordY1, int prmCoordX2, int prmCoordY2) { double[] mapScale = mapField.getSizeCell(); return Math.Round((Math.Sqrt(Math.Pow(prmCoordX2 * mapScale[0] prmCoordX1 * mapScale[0], 2) + Math.Pow(prmCoordY2 * mapScale[1] - prmCoordY1 * mapScale[1], 2))), 0); } } // Класс скважины нефтегазового месторождения class Well { public int numberWell; public int coordXWell; public int coordYWell; // пустой конструктор public Well() 122

{ } // конструктор public Well(int prmNumberWell, int prmCoordX, int prmCoordY) { numberWell = prmNumberWell; coordXWell = prmCoordX; coordYWell = prmCoordY; } } // Класс для работы с сеткой месторождения class GridField { private int sizeGridX; // количество ячеек по X private int sizeGridY; // количество ячеек по Y private double[,] gridPiezoConductivity; // значения коэффициента пьезопроводности в ячейках сетки // пустой конструктор public GridField() { } // Инициализация сетки месторождения по параметрам public void InitGridField(int prmSizeGridX, int prmSizeGridY, double prmValuePiezo) { sizeGridX = prmSizeGridX; sizeGridY = prmSizeGridY; gridPiezoConductivity = new double[sizeGridX, sizeGridY]; for (int x = 0; x < gridPiezoConductivity.GetLength(0); x++) for (int y = 0; y < gridPiezoConductivity.GetLength(1); y++) gridPiezoConductivity[x, y] = prmValuePiezo; } } Приложение Б.2. Программная реализация алгоритмов математической модели динамики эксплуатации нефтегазового месторождения. // // // // // // // public class Model { public DevField devField; public MapField mapField; // private int countProdWells; private int countInjWells; // история эксплуатации скважины public Classes.Parameters.ListModelHistoryRateWell listExplRateWell; история эксплуатации скважин // параметры скважин public Classes.Parameters.ListModelDistanceWells listParDistanceWells; расстояние между скважинами // параметрв среды public Classes.Parameters.Parameter parAint; параметр месторождения Aint public Classes.Parameters.Parameter parAout; параметр месторождения Aout // public Classes.Parameters.ListParDampingWell listParDampingWell; параметр затухания дебита (лямбда) public Classes.Parameters.ListParInterInjWell listParInterInjWell; параметр интерференции нагнетательной скважины public Classes.Parameters.ListParInterProdWell listParInterProdWell; параметр интерференции добывающей скважины 123

public Classes.Parameters.ListParSpeedFluidWells listParSpeedFluidInjWell; // параметр скорости движения флюида от нагнетательной скважины public Classes.Parameters.ListParSpeedFluidWells listParSpeedFluidProdWell; // параметр скорости движения флюида от добыающей скважины // public Classes.Parameters.ListModelRateWell listModelRateWell; // расчеты 1 public Classes.Parameters.listModelInjWell listModelInjWell; // влияние нагнетательных скважин public Classes.Parameters.listModelInjWell listModelProdWell; // влияние добывающих скважин // public Model() { devField = new DevField(); // Создание объекта месторождения //gridPool = new GridPool(); // Создание объекта сетки месторождения mapField = new MapField(); //iterationProcess = new IterationProcess(); CreateListParameters(); // инициализация структур для хранения параметров модели } // Расчет компоненты модели Q1 по номеру добывающей скважины и моменту времени private double CalcModelQ1(int prmNumberProdWell, int prmTime) { return Math.Exp( -(prmTime - 1) * listParDampingWell.GetParValue(prmNumberProdWell)) * listExplRateWell.GetParValueInit(prmNumberProdWell); } // Расчет компоненты модели Q2 по номеру добывающей скважины и моменту времени private double CalcModelQ2(int prmNumberProdWell, int prmTime) { double lQ2 = 0; for (int j = 0; j < devField.CountInjWells; j++) { int lNumberInjWell = listExplRateWell.listInjWells[j]; double lCoefLagTime = (listParDistanceWells.GetParValue(prmNumberProdWell, lNumberInjWell) / listParSpeedFluidInjWell.GetParValue(prmNumberProdWell, lNumberInjWell)); double lSum = 0; double lSumInjWater = 0; int lNumberTimeWithLag = 1; // счетчик для остлеживания момента времени с запаздыванием for (int t = 0; t < prmTime; t++) { int lLagTime = listExplRateWell.listTimeIntervals[t] (int)lCoefLagTime; // проверяем запаздывание if (lLagTime > 0) { lSum = lSum + Math.Abs(listExplRateWell.GetParValue(lNumberInjWell, lNumberTimeWithLag)) * 1; //* 1 - размер временного интервала lNumberTimeWithLag++; } } if (lSum > 0) lSumInjWater = listModelInjWell.GetParValueSumByInjWell(lNumberInjWell, prmTime); 124

double lSumInjWell = listParInterInjWell.GetParValue(prmNumberProdWell, lNumberInjWell) * (lSum lSumInjWater) * Math.Exp(-1 * parAint.GetParValue() * lCoefLagTime); listModelInjWell.AddParValue(lNumberInjWell, prmNumberProdWell, prmTime, lSumInjWell); lQ2 = lQ2 + lSumInjWell; } return lQ2; } // Расчет компоненты модели Q3 по номеру добывающей скважины и моменту времени private double CalcModelQ3(int prmNumberProdWell, int prmTime, int prmSizeTimeWindow) { double lQ3 = 0; for (int i = 0; i < devField.CountProdWells; i++) { double lSumProdWell = 0; int lNumberProdWellTo = listExplRateWell.listProdWells[i]; if (prmNumberProdWell != lNumberProdWellTo) { double lCoefLagTime = listParDistanceWells.GetParValue(prmNumberProdWell, lNumberProdWellTo) / listParSpeedFluidProdWell.GetParValue(prmNumberProdWell, lNumberProdWellTo); double lRateWellCurrent = 0; double lRateWellTo = 0; int lLagTime = (prmTime - 1) - (int)lCoefLagTime; // * if (lLagTime > 0) { if (lLagTime > prmSizeTimeWindow) // если рассматриваемое время превысило лимит { lRateWellCurrent = listModelRateWell.GetParValueSum(prmNumberProdWell, lLagTime); lRateWellTo = listModelRateWell.GetParValueSum(lNumberProdWellTo, lLagTime); lSumProdWell = listParInterProdWell.GetParValue(prmNumberProdWell, lNumberProdWellTo) * (lRateWellCurrent - lRateWellTo) * Math.Exp(-1 * parAout.GetParValue() * lCoefLagTime); } else // если рассматриваевое время входит в лимит { lRateWellCurrent = listExplRateWell.GetParValue(prmNumberProdWell, lLagTime); lRateWellTo = listExplRateWell.GetParValue(lNumberProdWellTo, lLagTime); lSumProdWell = listParInterProdWell.GetParValue(prmNumberProdWell, lNumberProdWellTo) * (lRateWellCurrent - lRateWellTo) * Math.Exp(-1 * parAout.GetParValue() * lCoefLagTime); } } listModelProdWell.AddParValue(prmNumberProdWell, lNumberProdWellTo, prmTime, lSumProdWell); lQ3 = lQ3 + lSumProdWell; } } return lQ3; } // расчет компонентов модели для эксперимента public void CalcModelforExp(int prmNumberProdWell, int prmNumberinjWell, int prmTime) { 125

listProdRateWellforExp.AddVoidPar(prmNumberProdWell, prmTime); listProdRateWellforExp.SetParValueQ1(prmNumberProdWell, prmTime, CalcModelQ1forExp(prmNumberProdWell, prmTime)); listProdRateWellforExp.SetParValueQ2(prmNumberProdWell, prmTime, CalcModelQ2forExp(prmNumberProdWell, prmNumberinjWell, prmTime)); listProdRateWellforExp.SetParValueQ3(prmNumberProdWell, prmTime, CalcModelQ3forExp(prmNumberProdWell, prmTime)); } public void CalcModelByProdWell(int prmNumberProdWell, int prmTime) { listModelRateWell.SetParValueQ1(prmNumberProdWell, prmTime, CalcModelQ1(prmNumberProdWell, prmTime)); listModelRateWell.SetParValueQ2(prmNumberProdWell, prmTime, CalcModelQ2(prmNumberProdWell, prmTime)); listModelRateWell.SetParValueQ3(prmNumberProdWell, prmTime, CalcModelQ3(prmNumberProdWell, prmTime, devField.CountTimeIntervals)); // listModelRateWell.CalcError(prmNumberProdWell, prmTime, listExplRateWell.GetParValue(prmNumberProdWell, prmTime)); } // Расчет целевой функции модели (полный) public double CalcObjectiveFunctionModel() { double lResult = 0; for (int t = 0; t < devField.CountTimeIntervals; t++) { int lTime = listExplRateWell.listTimeIntervals[t]; for (int i = 0; i < devField.CountProdWells; i++) { int lNumberProdWell = listExplRateWell.listProdWells[i]; lResult = lResult + Math.Abs(listExplRateWell.GetParValue(lNumberProdWell, lTime) listModelRateWell.GetParValueSum(lNumberProdWell, lTime)); } } return lResult; } // Расчет целевой функции модели на интервале времени public double CalcObjectiveFunctionModel(int prmTimeStart, int prmTimeEnd) { double lResult = 0; for (int t = prmTimeStart; t < prmTimeEnd; t++) { int lTime = listExplRateWell.listTimeIntervals[t]; for (int i = 0; i < devField.CountProdWells; i++) { int lNumberProdWell = listExplRateWell.listProdWells[i]; lResult = lResult + Math.Abs(listExplRateWell.GetParValue(lNumberProdWell, lTime) listModelRateWell.GetParValueSum(lNumberProdWell, lTime)); } } return lResult; } // Расчет целевой функции модели по скважине и моменту времени public double CalcObjectiveFunctionWellTime() { double lResult = 0; for (int t = 0; t < devField.CountTimeIntervals; t++) { int lTime = listExplRateWell.listTimeIntervals[t]; for (int i = 0; i < devField.CountProdWells; i++) 126

{ int lNumberProdWell = listExplRateWell.listProdWells[i]; lResult = lResult + Math.Abs(listExplRateWell.GetParValue(lNumberProdWell, lTime) listModelRateWell.GetParValueSum(lNumberProdWell, lTime)); } } return (lResult / (devField.CountProdWells * devField.CountTimeIntervals)); } // выгрузка параметров модели в линейный вектор public double[] GetVectorParamsHJ(int prmCountVectorParams) { double[] vectorParams = new double[prmCountVectorParams]; int k = 0; // счетчик по параметру vectorParams[k] = parAint.GetParValue(); // 1 k++; vectorParams[k] = parAout.GetParValue(); // 2 k++; for (int i = 0; i < devField.CountProdWells; i++) // 3 Лямбда { vectorParams[k] = listParDampingWell.GetParValueByIndex(i); k++; } for (int i = 0; i < devField.CountInjWells; i++) { vectorParams[k] = listParInterInjWell.GetParValueLossWaterByIndex(i); k++; } for (int i = 0; i < devField.CountInjWells * devField.CountProdWells; i++) // 5 Бетта { vectorParams[k] = listParInterInjWell.GetParValueByIndex(i); k++; } for (int i = 0; i < devField.CountProdWells * (devField.CountProdWells - 1); i++) // 6 Гамма { vectorParams[k] = listParInterProdWell.GetParValueByIndex(i); k++; } for (int i = 0; i < devField.CountInjWells * devField.CountProdWells; i++) // 6 V1 { vectorParams[k] = listParSpeedFluidInjWell.GetParValueByIndex(i); k++; } for (int i = 0; i < devField.CountProdWells * (devField.CountProdWells - 1); i++) // 6 V2 { vectorParams[k] = listParSpeedFluidProdWell.GetParValueByIndex(i); k++; } return vectorParams; } public double[] GetVectorParams(int prmCountVectorParams) { double[] vectorParams = new double[prmCountVectorParams]; int k = 0; // счетчик по параметру vectorParams[k] = parAint.GetParValue(); // 1 127

k++; vectorParams[k] = parAout.GetParValue(); // 2 k++; for (int i = 0; i < devField.CountProdWells; i++) // 3 Лямбда { vectorParams[k] = listParDampingWell.GetParValueByIndex(i); k++; } for (int i = 0; i < devField.CountInjWells * devField.CountProdWells; i++) // 4 Бетта { vectorParams[k] = listParInterInjWell.GetParValueByIndex(i); k++; } for (int i = 0; i < devField.CountInjWells * devField.CountProdWells; i++) // 5 V1 { vectorParams[k] = listParSpeedFluidInjWell.GetParValueByIndex(i); k++; } for (int i = 0; i < devField.CountProdWells * (devField.CountProdWells - 1); i++) // 6 Гамма { vectorParams[k] = listParInterProdWell.GetParValueByIndex(i); k++; } for (int i = 0; i < devField.CountProdWells * (devField.CountProdWells - 1); i++) // 7 V2 { vectorParams[k] = listParSpeedFluidProdWell.GetParValueByIndex(i); k++; } return vectorParams; } // загрузка параметров из линейного вектора в модель public void SetVectorParamsHJ(double[] prmVectorParams) { int k = 0; // счетчик по параметру parAint.SetParValue(prmVectorParams[k]); // 1 k++; parAout.SetParValue(prmVectorParams[k]); // 2 k++; for (int i = 0; i < devField.CountProdWells; i++) // 3 Лямбда { listParDampingWell.SetParValueByIndex(i, prmVectorParams[k]); k++; } for (int i = 0; i < devField.CountInjWells; i++) // 3 Losswater { k++; } for (int i = 0; i < devField.CountProdWells * devField.CountInjWells; i++) // 4 Бетта { listParInterInjWell.SetParValueByIndex(i, prmVectorParams[k]); k++; } for (int i = 0; i < devField.CountProdWells * (devField.CountProdWells - 1); i++) // 5 Гамма { listParInterProdWell.SetParValueByIndex(i, prmVectorParams[k]); 128

k++; } for (int i = 0; i < devField.CountProdWells * devField.CountInjWells; i++) // 6 V1 { listParSpeedFluidInjWell.SetParValueByIndex(i, prmVectorParams[k]); k++; } for (int i = 0; i < devField.CountProdWells * (devField.CountProdWells - 1); i++) // 6 V2 { listParSpeedFluidProdWell.SetParValueByIndex(i, prmVectorParams[k]); k++; } } public void SetVectorParams(double[] prmVectorParams) { int k = 0; // счетчик по параметру parAint.SetParValue(prmVectorParams[k]); // 1 k++; parAout.SetParValue(prmVectorParams[k]); // 2 k++; for (int i = 0; i < devField.CountProdWells; i++) // 3 Лямбда { listParDampingWell.SetParValueByIndex(i, prmVectorParams[k]); k++; } for (int i = 0; i < devField.CountProdWells * devField.CountInjWells; i++) // 4 Бетта { listParInterInjWell.SetParValueByIndex(i, prmVectorParams[k]); k++; } for (int i = 0; i < devField.CountProdWells * devField.CountInjWells; i++) // 5 V1 { listParSpeedFluidInjWell.SetParValueByIndex(i, prmVectorParams[k]); k++; } for (int i = 0; i < devField.CountProdWells * (devField.CountProdWells - 1); i++) // 6 Гамма { listParInterProdWell.SetParValueByIndex(i, prmVectorParams[k]); k++; } for (int i = 0; i < devField.CountProdWells * (devField.CountProdWells - 1); i++) // 7 V2 { listParSpeedFluidProdWell.SetParValueByIndex(i, prmVectorParams[k]); k++; } } Приложение Б.3. Программная реализация алгоритмов построения прогноза пространственного распределения коэффициента пьезопроводности проницаемого пласта. 129

class hydroTom { public int countPairWells; // количество пар скважин public int countSourceWells; // количество скважин источников сигнала public int countReceiverWells; // количество скважин приемников сигнала public List<intervalTimes> listIntervalTimes; // список временных интервалов public List<trackWells> listPairWells; public int[,] gridCountLineInPoint; // матрица содержит информацию про количество траекторий проходящих через точки сетки public double errorFull; private void CalcError() { errorFull = 0; for (int p = 0; p < countPairWells; p++) { errorFull += listIntervalTimes[p].errorTime; } errorFull = errorFull / countPairWells; } // конструктор public hydroTom(int x, int y) { listIntervalTimes = new List<intervalTimes>(); listPairWells = new List<trackWells>(); gridCountLineInPoint = new int[x, y]; } // Расчет времени и траектории для пары скважин (вариант 1) private void CalcTimeAndLengthPairWells(DevField.DevField prmDevField, Point prmReceiverWellPoint, int prmNumberPairWells, int x) { // Начало траектории (по Беллману ее окончание) DevField.Well sourceWell = prmDevField.GetWell(listIntervalTimes[prmNumberPairWells].numberSourceWell); // получение координаты скважины-источника сигнала Point pointSourceWell = new Point(sourceWell.coordXWell, sourceWell.coordYWell); Point pointCurrent = prmReceiverWellPoint; // рассматриваемая точка int countTracks = 0; // количество промежуточных траеткорий // номер точки в траектории int numbLevel = 1; // номер уровня траеткории // флаг выполнения процесса Point coordPoint = GetCoordLine(prmReceiverWellPoint, pointSourceWell); if (coordPoint.X != 0) { countTracks = prmDevField.mapField.sizeGridY; listPairWells[prmNumberPairWells].str = "Идем вдоль X"; } // количество промежуточных траекторий по оси X if (coordPoint.Y != 0) { countTracks = prmDevField.mapField.sizeGridX; listPairWells[prmNumberPairWells].str = "Идем вдоль Y"; } // количество промежуточных траекторий по оси Y // Запускам поиск оптимальной траектории по принципу Беллмана List<Point> lListPoint = GetListPoint(pointCurrent, coordPoint, countTracks); // получаем список точек для перехода /* * первый шаг по Беллману */ 130

X или Y) // строим переход из первой точки в каждую на соседней линии (по оси for (int i = 0; i < countTracks; i++) { listPairWells[prmNumberPairWells].AddTrack(pointCurrent); // создание траектории и добавление первой точки listPairWells[prmNumberPairWells].AddPointTrack(i, lListPoint[i], numbLevel); // добавление второй точки //double[] mas = CalcTimeAndLenght(prmMapField, pointCurrent, lListPoint[i], ); // расчет времени и длины double[] mas = CalcTimeAndLenght(prmDevField.mapField, pointCurrent, lListPoint[i], coordPoint); // расчет времени и длины listPairWells[prmNumberPairWells].AddTimeAndlength(i, mas[1], mas[0]); // сохранение значений времени и длины } pointCurrent = lListPoint[0]; numbLevel++; // Расчет количества переходов в каждой траекторий int countLevel = GetCountLevel(pointSourceWell, prmReceiverWellPoint, coordPoint); int tempCountLevel = countLevel; if (x != 0) { tempCountLevel = countLevel; countLevel = x; } while (numbLevel < countLevel) // пока не дошли до последней линии { track[] tempListTracksWell = new track[countTracks]; lListPoint = GetListPoint(pointCurrent, coordPoint, countTracks); // получаем список точек для перехода for (int i = 0; i < lListPoint.Count; i++) { int t = 0; // номер траектории для продолжения double optVariant = 9999999999999999; double optCurVariant = 0; // оптимальное время между рассматриваемой парой точек double[] masOpt = new double[2]; pointCurrent = lListPoint[i]; // точка из новой линии for (int numbT = 0; numbT < listPairWells[prmNumberPairWells].listTracks.Count; numbT++) { Point nextPoint = listPairWells[prmNumberPairWells].listTracks[numbT].GetLastPoint(); // точка для возможного перехода double[] mas = CalcTimeAndLenght(prmDevField.mapField, pointCurrent, nextPoint, coordPoint); // расчет времени и длины для двух точек double sumAll = listPairWells[prmNumberPairWells].GetSumTimeForPoints(numbT, 0, numbLevel); // получить время до текущей точки sumAll = Math.Round(sumAll + mas[1], 3); if (sumAll < optVariant) // новое время лучше текущего? { // да optVariant = sumAll; masOpt = mas; t = numbT; optCurVariant = mas[1]; } else { // нет 131

if (sumAll == optVariant) // новое и текущее время равно? { // да if (mas[1] < optCurVariant) // новый отрезок лучше текущего? { // да optVariant = sumAll; masOpt = mas; t = numbT; optCurVariant = mas[1]; //MessageBox.Show("Зашли"); } } } } tempListTracksWell[i] = new track(listPairWells[prmNumberPairWells].GetTrack(t)); tempListTracksWell[i].AddPointTrack(pointCurrent); tempListTracksWell[i].SetTimeAndLength(masOpt[1], masOpt[0]); } for (int i = 0; i < countTracks; i++) { listPairWells[prmNumberPairWells].listTracks[i] = tempListTracksWell[i]; } numbLevel++; } // переход в последнюю точку if (tempCountLevel == numbLevel) { pointCurrent = pointSourceWell; double optTime = 99999999; double curOptTime = 0; for (int numbT = 0; numbT < countTracks; numbT++) { Point nextPoint = listPairWells[prmNumberPairWells].listTracks[numbT].GetLastPoint(); // точка для возможного перехода listPairWells[prmNumberPairWells].AddPointTrack(numbT, pointCurrent, numbLevel); // добавление второй точки //double[] mas = CalcTimeAndLenght(prmMapField, pointCurrent, nextPoint); // расчет времени и длины double[] mas = CalcTimeAndLenght(prmDevField.mapField, pointCurrent, nextPoint, coordPoint); // расчет времени и длины для двух точек listPairWells[prmNumberPairWells].AddTimeAndlength(numbT, mas[1], mas[0]); // сохранение значений времени и длины double sum = listPairWells[prmNumberPairWells].GetSumTime(numbT); if (optTime > listPairWells[prmNumberPairWells].GetSumTime(numbT)) // выбор оптимальной траеткории { optTime = listPairWells[prmNumberPairWells].GetSumTime(numbT); listPairWells[prmNumberPairWells].timeOptTrack = optTime; listPairWells[prmNumberPairWells].lengthOptTrack = listPairWells[prmNumberPairWells].GetSumLength(numbT); listPairWells[prmNumberPairWells].numbOpt = numbT; curOptTime = mas[1]; 132

} else { if (optTime == listPairWells[prmNumberPairWells].GetSumTime(numbT)) // выбор оптимальной траеткории { if (mas[1] < curOptTime) { optTime = listPairWells[prmNumberPairWells].GetSumTime(numbT); listPairWells[prmNumberPairWells].timeOptTrack = optTime; listPairWells[prmNumberPairWells].lengthOptTrack = listPairWells[prmNumberPairWells].GetSumLength(numbT); listPairWells[prmNumberPairWells].numbOpt = numbT; curOptTime = mas[1]; } } } } // Время движения между точками // 0 - время, 1 - расстояние (вариант 1) private double[] CalcTimeAndLenght(DevField.MapField prmMapField, Point prmPoint1, Point prmPoint2, Point prmCoordPoint) { double[] mas = new double[2]; double[] masScale = prmMapField.getSizeCell(); // получаем размер ячейки в реальных единицах (м) double coefPiezo = CalcAvgPiezoconductivityPoints(prmMapField, prmPoint1, prmPoint2, prmCoordPoint); double lenghtTrack = CalcLineTrackBetweenPoints(prmPoint1, prmPoint2, masScale[0], masScale[1]); mas[0] = Math.Round(lenghtTrack, 2); // Траектория (длина) mas[1] = Math.Round((lenghtTrack / 3 * coefPiezo), 3); // Время return mas; } // Функция вовзращает среднее значение пьезопроводности области между выбранными ячейками (ячейки несмежные) область между ячейками выбирается в виде прямоугольника, который ограничивается этими ячейками private double CalcAvgPiezoconductivityPoints(DevField.MapField prmMapField, Point prmPointStart, Point prmPointEnd, Point prmCoordPoints) { double lAvgPiezo = 0; if ((Math.Abs(prmPointStart.X - prmPointEnd.X) <= 1) && (Math.Abs(prmPointStart.Y - prmPointEnd.Y) <= 1)) // являеются ли ячейки смежными { // да lAvgPiezo = CalcPiezoconductivityPoints(prmMapField.GetPiezoConductivityCellField(prmPointSt art.X, prmPointStart.Y), prmMapField.GetPiezoConductivityCellField(prmPointEnd.X, prmPointEnd.Y)); } else { // нет int type = 0; if (prmCoordPoints.X != 0) type = 1; // движение вдоль оси по X if (prmCoordPoints.Y != 0) type = 2; // движение вдоль оси по Y lAvgPiezo += prmMapField.GetPiezoConductivityCellField(prmPointStart.X, prmPointStart.Y); lAvgPiezo += prmMapField.GetPiezoConductivityCellField(prmPointEnd.X, prmPointEnd.Y); 133

switch (type) { case 1: // движение вдоль X (перебор ячеек по Y) { int count = Math.Abs(prmPointStart.Y prmPointEnd.Y); int k = 2; int step = 0; if (prmPointStart.Y > prmPointEnd.Y) step = 1; if (prmPointEnd.Y > prmPointStart.Y) step = -1; Point lPoint1 = new Point(prmPointStart.X, prmPointEnd.Y + step); Point lPoint2 = new Point(prmPointEnd.X, prmPointEnd.Y + step); for (int i = 1; i < count; i++) { lAvgPiezo += prmMapField.GetPiezoConductivityCellField(lPoint1.X, lPoint1.Y); lAvgPiezo += prmMapField.GetPiezoConductivityCellField(lPoint2.X, lPoint2.Y); k += 2; lPoint1.Y = lPoint1.Y + step; lPoint2.Y = lPoint2.Y + step; } lAvgPiezo = lAvgPiezo / k; break; } case 2: // движение вдоль Y (перебор ячеек по X) { int step = 0; if (prmPointStart.X < prmPointEnd.X) step = -1; else step = 1; Point lPoint1 = new Point(prmPointEnd.X + step, prmPointStart.Y); Point lPoint2 = new Point(prmPointEnd.X + step, prmPointEnd.Y); int count = Math.Abs(prmPointStart.X prmPointEnd.X); int k = 0; for (int i = 1; i < count; i++) { lAvgPiezo += prmMapField.GetPiezoConductivityCellField(lPoint1.X, lPoint1.Y); lAvgPiezo += prmMapField.GetPiezoConductivityCellField(lPoint2.X, lPoint2.Y); k += 2; lPoint1.X = lPoint1.X + step; lPoint2.X = lPoint2.X + step; } lAvgPiezo = lAvgPiezo / k; break; } } } return lAvgPiezo; } // Функция рассчитывает и возвращает коэффициент пьезопроводности между точками сетки private double CalcPiezoconductivityPoints(double prmPiezoPointStart, double prmPiezoPointEnd) { return (prmPiezoPointStart + prmPiezoPointEnd) / 2; } 134

сетки // Функция рассчитывает и возвращает длину траектории между точками private double CalcLineTrackBetweenPoints(Point prmPointStart, Point prmPointEnd, double prmSizeX, double prmSizeY) { return Math.Round(Math.Sqrt(Math.Pow(prmPointEnd.X * prmSizeX prmPointStart.X * prmSizeX, 2) + Math.Pow(prmPointEnd.Y * prmSizeY prmPointStart.Y * prmSizeY, 2)), 2); } // Функция расчитывает и возвращает длину участка между точками (отдельно) public double CalcLength(DevField.MapField prmMapField, Point prmPoint1, Point prmPoint2) { double[] masScale = prmMapField.getSizeCell(); return CalcLineTrackBetweenPoints(prmPoint1, prmPoint2, masScale[0], masScale[1]); } // Функция расчитывает количество линий в точке public void CalcLineInPoint() { for (int i = 0; i < countPairWells; i++) { int numbOpt = listPairWells[i].numbOpt; for (int n = 0; n < listPairWells[i].listTracks[numbOpt].listPointsTrack.Count; n++) { Point optPoint = listPairWells[i].listTracks[numbOpt].listPointsTrack[n].point; gridCountLineInPoint[optPoint.X - 1, optPoint.Y - 1] += 1; } } } } Приложение Б.4. Программная реализация алгоритмов построения нечеткой модели фильтрационного сопротивления проницаемого пласта по результатам работы гидродинамической томографии class fuzzyElements { public int sizeX; public int sizeY; List<Point> listPoints; double[,] fieldScattering; public double[,] fieldInterpolation; double[,] functionMemberShip; double parameterEffective; public int[,] temp; // Функция интерполяции public void Interpolation() { temp = new int[sizeX, sizeY]; fieldInterpolation = new double[sizeX, sizeY]; double interpolationVal = 0; for (int x = 0; x < sizeX; x++) 135

{ for (int y = 0; y < sizeY; y++) { fieldInterpolation[x, y] = fieldScattering[x, y]; } } int countIt = 0; while (CheckNullCell(fieldInterpolation) != 0) { List<el> list = new List<el>(); for (int x = 0; x < sizeX; x++) { for (int y = 0; y < sizeY; y++) { list.Add(new el(x, y, fieldInterpolation[x, y])); //fieldInterpolation[x, y] = fieldScattering[x, y]; } } countIt++; temp = new int[sizeX, sizeY]; //MessageBox.Show("Количество пустых яччеек: " + CheckNullCell(fieldInterpolation).ToString()); for (int x = 0; x < sizeX; x++) { for (int y = 0; y < sizeY; y++) { interpolationVal = fieldInterpolation[x, y]; if (interpolationVal == 0) { for (int x2 = x - 1; x2 <= x + 1; x2++) { for (int y2 = y - 1; y2 <= y + 1; y2++) { if (((x2 > 0) && (x2 < sizeX)) && ((y2 > 0) && (y2 < sizeY))) { if (list.Find(par => par.x == x2 && par.y == y2).val > 0) temp[x, y] += 1; } } } } } } int max = FindMax(temp); //MessageBox.Show("Максимум: " + max.ToString()); for (int x = 0; x < sizeX; x++) { for (int y = 0; y < sizeY; y++) { double valInter = 0; int count = 0; if ((temp[x, y] == max) && (max != 0)) { for (int x2 = x - 1; x2 <= x + 1; x2++) { for (int y2 = y - 1; y2 <= y + 1; y2++) { if (((x2 > 0) && (x2 < sizeX)) && ((y2 > 0) && (y2 < sizeY))) { if (fieldInterpolation[x2, y2] != 0) { 136

valInter += fieldInterpolation[x2, y2]; count++; } } } } if (count != 0) { valInter = valInter / count; fieldInterpolation[x, y] = valInter; } } } } if (countIt == 18) break; } } // Функция возвращает количество нулевых ячеек private int CheckNullCell(double[,] prmGrid) { int countNullCell = 0; for (int x = 0; x < sizeX; x++) { for (int y = 0; y < sizeY; y++) { if (prmGrid[x, y] == 0) countNullCell++; } } return countNullCell; } // Функция возращает значения максимума по сетке private int FindMax(int[,] prmGrid) { int max = 0; for (int x = 0; x < sizeX; x++) for (int y = 0; y < sizeY; y++) (temp[x, y] > max) max = temp[x, y]; return max; } if public double ParameterEffective {get { return parameterEffective; } set { parameterEffective = value; }} // // конструктор public fuzzyElements(int prmSizeX, int prmSizeY) { sizeX = prmSizeX; sizeY = prmSizeY; listPoints = new List<Point>(); } // Функция загрузки списка точек public void FillListPoints(HDT.hydroTom prmHT) { for (int i = 0; i < prmHT.countPairWells; i++) { for (int n = 0; n < prmHT.listPairWells[i].listTracks[prmHT.listPairWells[i].numbOpt].listPointsTrac k.Count - 1; n++) { if (!listPoints.Contains(prmHT.listPairWells[i].listTracks[prmHT.listPairWells[i].n umbOpt].listPointsTrack[n].point)) 137

listPoints.Add(prmHT.listPairWells[i].listTracks[prmHT.listPairWells[i].numbOpt] .listPointsTrack[n].point); } } } //Функция загрузки данных по точках из модели public void FillScatteringField(double[,] prmGridPiezo) { fieldScattering = new double[sizeX, sizeY]; for (int x = 0; x < sizeX; x++) { for (int y = 0; y < sizeY; y++) { if (listPoints.Contains(new Point(x, y))) fieldScattering[x - 1, y - 1] = prmGridPiezo[x, y]; else fieldScattering[x, y] = 0; } } } public double[,] GetFieldScattering() { return fieldScattering; } public double[,] GetFieldMemberShip() { return functionMemberShip; } // расчет функции принадлежности с суммой (без учета матрицы траекторий) public void CalcSumFunctionMemberShip(double prmParamterEffective) { parameterEffective = prmParamterEffective; functionMemberShip = new double[sizeX, sizeY]; for (int i = 1; i <= sizeX; i++) { for (int j = 1; j <= sizeY; j++) { double sumFuncMShip = 0; double currentCoefPiezo = fieldScattering[i - 1, j - 1]; // значение пьезопроводности для которого рассчитываем функцию принадлежности for (int i2 = 1; i2 <= sizeX; i2++) { for (int j2 = 1; j2 <= sizeY; j2++) { double valFuncMShip = 0; double localPiezo = fieldScattering[i2 - 1, j2 - 1]; // формула оригинал if (localPiezo == 0) valFuncMShip = 0; else { double r; if ((i == i2) && (j == j2)) r = 1; else r = Math.Pow(i - i2, 2) + Math.Pow(j - j2, 2); //r = Math.Sqrt(Math.Pow(i - i2, 2) + Math.Pow(j - j2, 2)); // считаем расстояние между рассматриваемыми точками //valFuncMShip = 1 / (parameterEffective * Math.Sqrt(2 * Math.PI) * r) * Math.Exp(-1 * (Math.Pow(currentCoefPiezo localPiezo, 2)) / (2 * Math.Pow(parameterEffective, 2))); valFuncMShip = 1 / (parameterEffective * Math.Sqrt(2 * Math.PI)) * Math.Exp(-1 * (Math.Pow(currentCoefPiezo - localPiezo, 2)) / (2 * Math.Pow(parameterEffective, 2))) * Math.Exp(-1 * r / (2 * Math.Pow(parameterEffective, 2))); // первоначальная формула без корня в расстоянии Math.Exp(-1 * (Math.Pow(i - i2, 2) + Math.Pow(j - j2, 2)) 138

} sumFuncMShip += valFuncMShip; } } functionMemberShip[i - 1, j - 1] = sumFuncMShip; // значение функции принадлежности } } NormalFunctionMShip(); } // расчет функции принадлежности с верхней гранью (без учета матрицы траекторий) public void CalcMaxFunctionMemberShip(double prmParamterEffective) { parameterEffective = prmParamterEffective; functionMemberShip = new double[sizeX, sizeY]; for (int i = 1; i <= sizeX; i++) { for (int j = 1; j <= sizeY; j++) { //double maxFuncMShip = 0; double currentCoefPiezo = fieldScattering[i - 1, j - 1]; // значение пьезопроводности для которого рассчитываем функцию принадлежности double[,] pointMemShip = new double[sizeX, sizeY]; for (int i2 = 1; i2 <= sizeX; i2++) { for (int j2 = 1; j2 <= sizeY; j2++) { double valFuncMShip = 0; double localPiezo = fieldScattering[i2 - 1, j2 - 1]; // формула оригинал if (localPiezo == 0) valFuncMShip = 0; else { double r; if ((i == i2) && (j == j2)) r = 1; else r = Math.Sqrt(Math.Pow(i - i2, 2) + Math.Pow(j - j2, 2)); // считаем расстояние между точками valFuncMShip = 1 / (parameterEffective * Math.Sqrt(2 * Math.PI) * r) * Math.Exp(-1 * (Math.Pow(currentCoefPiezo localPiezo, 2)) / (2 * Math.Pow(parameterEffective, 2)));// * Math.Exp(-1 * / (2 * Math.Pow(parameterEffective, 2))); } //if (maxFuncMShip < valFuncMShip) maxFuncMShip = valFuncMShip; pointMemShip[i2 - 1, j2 - 1] = valFuncMShip; } } functionMemberShip[i - 1, j - 1] = FindMax(pointMemShip); //maxFuncMShip; // значение функции принадлежности } } //NormalFunctionMShip(); } // расчет функции принадлежности с нижней гранью* под большим вопросом public void CalcMinFunctionMemberShip(double prmParamterEffective) { parameterEffective = prmParamterEffective; functionMemberShip = new double[sizeX, sizeY]; for (int i = 1; i <= sizeX; i++) { for (int j = 1; j <= sizeY; j++) { 139

double minFuncMShip = 0; double currentCoefPiezo = fieldScattering[i - 1, j - 1]; // значение пьезопроводности для которого рассчитываем функцию принадлежности for (int i2 = 1; i2 <= sizeX; i2++) { for (int j2 = 1; j2 <= sizeY; j2++) { double valFuncMShip = 0; double localPiezo = fieldScattering[i2 - 1, j2 - 1]; if (localPiezo == 0) valFuncMShip = 0; else valFuncMShip = (1 / parameterEffective * Math.Sqrt(2 * Math.PI)) * Math.Exp(-1 * (Math.Pow(currentCoefPiezo - localPiezo, 2)) / (2 * Math.Pow(parameterEffective, 2))) * Math.Exp(-1 * (Math.Pow(i - i2, 2) + Math.Pow(j - j2, 2)) / (2 * Math.Pow(parameterEffective, 2))); if (minFuncMShip > valFuncMShip) minFuncMShip = valFuncMShip; } } functionMemberShip[i - 1, j - 1] = minFuncMShip; // значение функции принадлежности } } NormalFunctionMShip(); } private double FindMin(double[,] prmField) { double min = 1;// prmField[0, 0]; for (int x = 0; x < prmField.GetLength(0); x++) for (int y = 0; y < prmField.GetLength(1); y++) if (prmField[x, y] != 0) if (prmField[x, y] < min) min = prmField[x, y]; return min; } private double FindMax(double[,] prmField) { double max = 0; // prmField[0, 0]; for (int x = 0; x < prmField.GetLength(0); x++) for (int y = 0; y < prmField.GetLength(1); y++) if (max < prmField[x, y]) max = prmField[x, y]; return max; } private double FindSum(double[,] prmField) { double sum = 0; // prmField[0, 0]; for (int x = 0; x < prmField.GetLength(0); x++) for (int y = 0; y < prmField.GetLength(1); y++) sum += prmField[x, y]; return sum; } private double[,] NormalizationField(double[,] prmField, double prmNormValue) { for (int x = 0; x < prmField.GetLength(0); x++) for (int y = 0; y < prmField.GetLength(1); y++) prmField[x, y] = prmField[x, y] / prmNormValue; return prmField; } // public void CalcFunctionMemberShip2(double prmParamterEffective) { parameterEffective = prmParamterEffective; functionMemberShip = new double[sizeX, sizeY]; for (int i = 1; i <= sizeX; i++) { for (int j = 1; j <= sizeY; j++) { double sumFuncMShip = 0; 140

double currentCoefPiezo = fieldScattering[i - 1, j - 1]; // значение пьезопроводности для которого рассчитываем функцию принадлежности for (int i2 = 1; i2 <= sizeX; i2++) { for (int j2 = 1; j2 <= sizeY; j2++) { double valFuncMShip = 0; //if ((i == i2)&(j == j2)) valFuncMShip = 0; //else { valFuncMShip = fieldScattering[i2 - 1, j2 - 1] * Math.Exp(-1* (Math.Pow(i - i2, 2) + Math.Pow(j - j2, 2)) / (2 * Math.Pow(parameterEffective, 2))); } sumFuncMShip += valFuncMShip; } } functionMemberShip[i - 1, j - 1] = sumFuncMShip; // значение функции принадлежности } } NormalFunctionMShip(); } // нормировка функции принадлежности private void NormalFunctionMShip() { double max = FindMax(); if (max != 0) for (int i = 0; i < sizeX; i++) for (int j = 0; j < sizeY; j++) functionMemberShip[i, j] = functionMemberShip[i, j] / max; } // поиск максимума private double FindMax() { double max = functionMemberShip[0, 0]; for (int x = 0; x < sizeX; x++) for (int y = 0; y < sizeY; y++) if (max < functionMemberShip[x, y]) max = functionMemberShip[x, y]; return max; 141

Рецензии:

Авторизуйтесь, чтобы добавить рецензию

- у работы пока нет рецензий -

Отзывы:

Авторизуйтесь, чтобы оставить отзыв