Математич еское модел ирование и ч исл енное прогнозирование характеристик природных динамич еских систем
На правах рукописи
Середа Альгирдас-Владимир Игнатьевич МАТЕМАТИЧ ЕСКОЕ МОДЕЛ ИРОВАНИЕ И Ч ИСЛ ЕННОЕ ПРОГНОЗИРОВАНИЕ ХАРАКТЕРИСТИК ПРИРОДНЫХ ДИНАМИЧ ЕСКИХ СИСТЕМ Специальность: 05.13.18 – Математическое моделирование, численные методы и комплексы программ
АВТОРЕФЕРАТ
диссертации на соискание ученой степени доктора технических наук
Санкт-Петербург – 2009 2
Работа выполнена в Мурманском государственном техническом университете
Официальные оппоненты:
доктор технических наук, профессор Бычков Ю.А.
доктор технических наук, профессор Истомин Е.П..
доктор физико-математических наук, профессор Клячкин В.И.
Ведущая организация – Санкт-Петербургский государственный морской технический университет
Защита диссертации состоится «_» 2009 года в _ часов на заседании совета по защите докторских и кандидатских диссертаций Д 212.238.01 Санкт Петербургского государственного электротехнического университета «ЛЭТИ» им.
В.И. Ульянова (Ленина) по адресу: 197376, Санкт-Петербург, ул. Проф. Попова, 5.
С диссертацией можно ознакомиться в библиотеке университета
Автореферат разослан «_» _ 200_ г.
Ученый секретарь совета
Общая характеристика работы
Актуальность работы. Под природной динамической системой (t) в данной работе понимается сплошная, пространственно-неоднородная по своим свойствам и составу, структурно определенная, в общем случае многокомпонентная материальная среда. Система (t) формируется во времени и заполняет собой некоторую ограниченную пространственную область (t), t[0,T], где t=0 – момент зарождения системы, t=T – настоящее время.
Формирование системы (t) включает в себя в общем случае процессы добавления в нее извне новых или исключения из нее части уже входивших в систему структурных элементов, а также изменение их свойств, вследствие протекающих в системе внутренних процессов различной природы (механических, физических, химических и т.д.). В общем случае структура системы и пространственные границы области (t) изменяются во времени.
Состояние системы (t) в каждый момент времени t[0,T] определяется ее структурой и физическими полями - пространственными распределениями в области (t) количественных значений физических характеристик, определяющих свойства образующей систему среды и присущих ей внутренних процессов. В общем случае текущее состояние системы является результатом ее эволюционного развития в предшествующий промежуток времени. Предполагается, что в силу объективных причин, структура системы (t) и пространственные распределения ее физических характеристик потенциально доступны для прямых или косвенных измерений лишь в настоящий момент времени t=T.
Изучение природных динамических систем с помощью непосредственных наблюдений и измерений, проведения натурных экспериментов и т.п., как правило, либо принципиально невозможно, либо сопряжено с большими временными и материальными затратами. В этой связи важными, а иногда и единственно возможными, инструментами исследования таких систем являются математическое моделирование и численный анализ. Вместе с тем, несмотря на очевидные успехи, использование математического аппарата и вычислительной техники в ряде случаев оказывается недостаточно эффективным с точки зрения прикладных целей исследования. Причиной этому служат многообразие и сложность протекающих в природных системах взаимосвязанных процессов, недостаток фактических априорных данных об условиях существования и свойствах исследуемых систем. В результате попытки точного описания приводят к чрезвычайно сложным для анализа математическим моделям, а недостаток данных не позволяет осуществлять адекватные реальным процессам вычислительные эксперименты. Упрощенные математические описания зачастую непригодны из-за большой погрешности, обусловленной игнорированием многих сопутствующих моделируемым процессам факторов и в силу объективно свойственной этим описаниям различного рода усредненности.
Возникает необходимость разработки методологии практически эффективного численного исследования природных динамических систем, с одной стороны позволяющего в достаточной степени учитывать их специфические особенности и свойства, а с другой стороны ориентированного на использование достаточно простых математических моделей и методов. Разработка такой методологии должна осуществляться с позиций системного подхода и носить целенаправленный характер, а ее использование должно приводить к рациональным вычислительным схемам, сбалансированным по точности с точностью и полнотой имеющихся данных о состояниях и свойствах изучаемых систем и процессов. Важным элементом результирующей компьютерной технологии должен быть проблемно ориентированный человеко-машинный интерфейс, обеспечивающий специалисту в соответствующей предметной области возможность анализа и корректировки хода вычислительного процесса на всех его этапах.
Одной из задач при исследовании природных динамических систем является задача определения значений, которые принимает в настоящий момент времени та или иная (целевая) физическая характеристика системы (t), в некоторой заданной подобласти V0 области (t), недоступной для прямых измерений. Такую задачу будем называть задачей прогнозирования характеристик природной динамической системы (t). Часто подобные задачи решаются посредством подходящей интерполяции или экстраполяции имеющихся данных о значениях целевой характеристики в требуемую пространственную область. Однако в ряде случаев такой подхо д к решению задачи не приемлем. В частности, это может иметь место в силу зависимости пространственного распределения целевой характеристики не только от конкретного состояния системы, но и от истории ее предшествовавшего развития. В этом случае задача должна решаться в контексте эволюционного развития системы от момента ее зарождения до настоящего времени, что является существенным отличием от традиционных задач обработки и интерпретации данных наблюдений. Разработка методологии практически эффективного решения задачи численного прогнозирования характеристик природных динамических систем, безусловно, является актуальной проблемой.
Естественным и практически важным примером природной динамической системы, в смысле приведенного выше описания, является геофлюидодинамическая система (ГФС), формирование которой происходит в осадочных бассейнах земли.
Многокомпонентность материальной среды ГФС предполагает, в том числе, ее многофазность, то есть, наличие в каждом элементарном объеме этой среды, как твердой, так и более подвижной материальной субстанции – флюида (жидкости или газа). Одной из важнейших характеристик ГФС является пространственное распределение геофлюидальных давлений 1 (ГФД) в осадочном бассейне. Знание ГФД в заданных пространственных областях имеет большое практическое значение, в частности, при бурении разведочных или промысловых скважин. Актуальность разработки методологии практически эффективного численного прогнозирования значений ГФД не вызывает сомнений. В данной работе эта методология разрабатывается применительно к задачам региональных разведочных работ в контексте предложенной общей методологии численного прогнозирования характеристик природных динамических систем.
В целом работа относится к области прикладного математического моделирования и численных методов в естественнонаучных задачах.
Целью диссертационной работы является разработка общей схемы методологии численного прогнозирования характеристик природных динамических систем и разработка на ее основе методологии численного прогнозирования ГФД.
Достижение указанной цели предполагает решение следующих задач.
давлений, воздействующих на фл юид в осадочных толщах земли.
1. Формулировка основных (концептуальных) принципов и построение на их основе общей схемы методологии численного прогнозирования характеристик природных динамических систем, включающей в себя:
- выбор, построение и анализ совокупности структурных и математических моделей исследуемой системы и протекающих в ней процессов;
- постановку, рациональное упрощение и численное решение соответствующих прямых и обратных задач;
- формирование и уточнение прогноза значений целевых характеристик исследуемой системы.
2. Постановка обратной задачи и разработка численных методов ее решения.
3. Разработка алгоритмических процедур рационального упрощения используемых моделей.
4. Разработка алгоритмической процедуры формирования прогноза значений целевой характеристики в заданной пространственной области.
5. Разработка алгоритмической процедуры уточнения прогноза.
6. Разработка, на основе общей схемы методологии численного прогнозирования, методологии численного прогнозирования ГФД в осадочных бассейнах земли, предполагающая, в дополнение к названным выше, решение следующих задач:
6.1. выбор, построение и анализ эффективной структурной (геологической) модели ГФС и эффективной математической модели флюидодинамического (ФД) процесса, декомпозиция математической модели;
6.2. постановка прямых задач;
6.3. разработка численных методов решения прямых задач;
6.4. конкретизация с учетом специфики ГФС алгоритмических процедур рационального упрощения используемых моделей, прогнозирования и уточнения прогноза значений ГФД в заданной области.
Методы исследования. В работе используется аппарат теории численных методов линейной алгебры, методов оптимизации и методов интерполирования данных, разностных методов решения краевых задач математической физики, теории регуляризирующих методов решения обратных задач.
Научная новизна работы заключается в следующем:
1. Сформулированы основные (концептуальные) принципы построения и разработана общая схема методологии численного прогнозирования характеристик природных динамических систем, основанная на системном подходе и эффективном моделировании в контексте эволюционного развития системы. Методология предусматривает оперативную оценку и корректировку хода вычислительного процесса на основе проблемно ориентированного человеко-машинного интерфейса.
В рамках общей схемы методологии численного прогнозирования предложены:
- численный метод решения обратных задач, основанный на концепции построения так называемых -квазирешений, позволяющий осуществлять параметрическую идентификацию (калибровку) соответствующих прямых задач;
- алгоритмическая процедура согласованной (региональной) калибровки прямых задач;
- алгоритмические процедуры рационального снижения размерности пространства модельных параметров для используемых эффективных моделей;
- алгоритмическая процедура формирования прогноза для заданной области;
- алгоритмическая процедура уточнения прогноза.
2. Разработана методология численного прогнозирования ГФД, представляющая собой конкретизацию общей схемы методологии численного прогнозирования применительно к ГФС. В рамках этой методологии:
- предложен комплекс эффективных математических моделей флюидодинамического процесса. В качестве основной использована математическая модель процесса фильтрации однофазного несжимаемого флюида в осадочном бассейне земли. Соответствующая прямая задача представляет собой краевую задачу для линейного параболического уравнения математической физики с переменными коэффициентами, исследуемую в геологическом масштабе времени (это позволяет отнести ее к задачам так называемого бассейнового моделирования) в трехмерной пространственной области с изменяющейся (по оси OZ) границей. В результате декомпозиции исходной математической модели по пространственным переменным предложены эффективные математические модели в одномерной («вертикальная» модель) и дву хмерной («латеральная» модель) пространственных областях.
Сформулированы соответствующие прямые задачи;
- предложены численные методы решения вертикальной и латеральной прямых задач;
- предложен метод совместного решения латеральной и вертикальной задач.
Практическая ценность работы заключается в том, что ней разработана методология численного прогнозирования ГФД в осадочных бассейнах земли.
Предлагаемый в работе общий (концептуальный) подход к построению методологии численного решения важной естественнонаучной задачи прогнозирования характеристик природных динамических систем может быть практически полезен при исследовании различных природных динамических систем.
Разработанные в рамках этой методологии алгоритмические решения носят общий характер и могут быть полезны при решении и других прикладных задач.
Достоверность научных и практических результатов. Достоверность научных положений и рекомендаций, приведенных в диссертации, подтверждается их достаточно строгим и аргументированным обсуждением в работе. Достоверность алгоритмических решений и численных методов практически подтверждается результатами их использования в расчетах с синтетическими и реальными данными.
Практическая эффективность методологии численного прогнозирования ГФД подтверждается результатами ее использования при обработке реальных данных в различных осадочных бассейнах.
Реализация результатов работы. На начальной стадии в 1994-1997 годах работа велась в НИИ МОРГЕОфизика (г.Мурманск) при поддержке Norsk Hydro (Берген, Норвегия). В 1994 году были разработаны концептуальные основы и первоначальные алгоритмические решения, положенные в основу компьютерной технологии PANDA 1 ® прогнозирования ГФД в осадочных бассейнах земли, разрабатывавшейся в 1995-1997 годах при поддержке Norsk Hydro и фирмы WST (Берген).
Pressure Analysis for Drilling Applications Well Service Technology.
В 1997–1999 годах в НИИ МОРГЕОфизика (г.Мурманск) выполнялся научно исследовательский проект в области бассейнового моделирования и прогнозирования ГФД.
В 1998-1999 годах работы проводились в рамках европейского научно исследовательского проекта JIP ODS 1 Project.
В 1999-2001 годах при поддержке НИИ МОРГЕОфизика Сервис.Ком (г.Мурманск) разрабатывалась вторая версия компьютерной технологии - “ПАНДА 2000©”.
В 2002—2005 годах исследования были поддержаны трехгодичным грантом АФГИР2 на основании трехлетнего договора между Мурманским государственным техническим университетом (МГТУ) и фирмой «Шлюмберже3 ».
Основные результаты работы (концепция, методология, модели и алгоритмы) реализованы в настоящее время в промышленном программном пакете «Панда 2000©», владельцем которого является НИИ МОРГ ЕОфизика Сервис.Ком, а также в исследовательских программных разработках, осуществлявшихся по результатам исследований в рамках гранта АФГИР. К ним относятся, в частности, пилотные версии программ согласованной (региональной) калибровки, декомпозиционных алгоритмов моделирования ФД процессов и ряд других.
По результатам работы получен патент Российской Федерации на изобретение №2321064 «Способ построения обратимой трехмерной гидродинамической модели земли, калибруемой в реальном масштабе времени в процессе бурения».
Основные элементы методологии численного прогнозирования ГФД, включая уточнение прогноза в процессе бурения, были успешно практически апробированы на реальных данных, полученных на территориях нефтегазовых месторождений России и других стран. В том числе:
В 1995-1998 годах. Викинг и Центральный грабен Северного моря.
1997 год. Медынская площадь Тимано-Печорской нефтегазовой провинции (НПГ).
2001-2003 годы. Площадь Ямбургского газового месторождения.
Научные положения, выносимые на защиту:
1. Концепция построения и общая схема методологии численного прогнозирования характеристик природных динамических систем, базирующиеся на принципах эволюционного развития, эффективного моделирования, минимальной сложности и пространственной локализации.
2. Численный метод решения обратных задач, основанный на концепции построения так называемых -квазирешений.
3. Алгоритмическая процедура согласованной (региональной) калибровки прямых задач.
4. Алгоритмические процедуры определения рационального количества структурных элементов и модельных параметров в используемых при исследовании системы эффективных структурных и математических моделях.
JIP ODS Project - The University of Liverpool, Fault Analysis Group, United Kingdom;
The University of Reading, Postgraduate Research Institute of Sedimentology, United Kingdom;
Osservatorio Geofisico Sperimentale, Italy;
Norsk Hydro, E & P Research Centre, Norway.
Американский фонд гражданских исследований и развития (CRDF) Schlumberger Research & Development Inc.
5. Алгоритмические процедуры формирования и уточнения прогноза значений целевой характеристики системы для заданной пространственной области.
6. Методология численного прогнозирования ГФД в осадочных бассейнах земли. В том числе:
- выбор эффективной математической модели флюидодинамического процесса и ее декомпозиция;
-постановка и численные методы решения прямых задач.
Апробация работы. Результаты представленных в работе исследований многократно докладывались на российских и международных научных и научно практических конференциях, в том числе на: международной научно-технической конференции (МНТК) "HPHT - Prediction", Workshop, Stavanger, Norway, 1994;
МНТК "57-th EA GE conference", Glasgow, Scotland, 1995;
МНТК "SEG, EA GO and EA GE international geophysical conference", St.Petersburg, 1995;
МНТК “Compaction and Overpressure Current Research”, Workshop, Institute Francais du Petrole, Paris, 1996;
МНТК "NORTHERN UNIVERSITIES", Murmansk, 1997;
МНТК “New methods and technologies in petroleum geology, drilling and reservoir engineering”, Workshop, Krakow, Poland, 1997;
МНТК «Pressure regimes in sedimentary basins and their prediction», American Association of Drilling Engineers (AADE) Forum, Houston, Texas, USA, 1998;
МНТК, посвященной 50-летию МГТУ, Мурманск, 2000;
всероссийской научно технической конференции (ВНТК) "Наука и образование - 2002", Мурманск, 2002;
ВНТК "Наука и образование - 2003", Мурманск, 2003;
МНТК "Наука и образование 2005", Мурманск, 2005;
МНТК "Наука и образование - 2006", Мурманск, 2006.
Публикации. По теме диссертации опубликована 21 работа, из них – 19 статей (12 статей опубликованы в рецензируемых научных журналах и изданиях, рекомендованных ВАК Минобрнауки России), одно рекламно-техническое описание (РТО) и один патент РФ на изобретение.
Структура и объем диссертации Диссертация состоит из введения, 6-ти глав, заключения, списка литературы, включающего 149 наименований. Основная часть работы изложена на 252 машинописных страницах. Работа содержит 39 рисунков.
Содержание работы Во введении содержится общая характеристика исследуемой проблемы, сформулирована цель работы, дано обоснование ее актуальности, приведено краткое изложение выполненных исследований и их результатов.
В первой главе вводятся необходимые понятия и определения, связанные с природной динамической системой. Приводится формальная постановка задачи прогнозирования характеристик природной динамической системы и осуществляется построение общей схемы методологии ее численного решения.
Вопросы теоретического и численного моделирования сложных динамических систем, в том числе вопросы решения некорректных задач отражены в работах Льюнга Л., Лебедева А.Н., Тихонова А.Н., Лаврентьева М.М., Алифанова О.М., Бакушинского А.Б., Гончарского А.В. и многих других о течественных и зарубежных авторов. Предлагаемая в данной работе методология разрабатывалась исходя из необходимости обеспечения возможности решения задачи прогнозирования характеристик природных динамических систем в реальном масштабе времени, что обеспечивается в результате построения и использования эффективных моделей исследуемых процессов, их целесообразной декомпозиции и рационального упрощения в соответствии с рядом дополнительных концептуальных положений (принципов), положенных в основу разработки общей схемы методологии. Важными составляющими методологии являются процедуры параметрической идентификации эффективных моделей исследуемых процессов, предполагающие постановку и решение обратных задач. В силу некорректности последних для их решения предлагается использовать регуляризирующий подход, ориентированный на построение так называемых -квазирешений. Сложность однозначной интерпретации имеющихся данных, оценки результатов расчетов на различных этапах вычислительного процесса предполагает включение в общую схему методологии проблемно ориентированного человеко-машинного интерфейса.
Рассмотрим идеализированное описание природной динамической системы.
Пусть (t) – природная динамическая система заполняет пространственную область (t), которая в любой момент времени t[0,T] может быть представлена как объединение пространственно упорядоченной совокупности конечного числа ограниченных, пространственно протяженных, односвязных и не имеющих общих точек подобластей i (t)(t), i=1,2,…,n (t), n (t)1:
n (t) (t) = i (t) ;
i (t)j (t) = ij, t[0,T].
(1) i= Каждая из подобластей заполнена физически однородной1 материальной средой. Будем называть эти подобласти слоями. Пространственная упорядоченность слоев i (t), i=1,2,…,n (t) определяет структуру исследуемой системы, а их количество и расположение в пространстве в общем случае зависят от времени.
Систематизированное описание пространственно упорядоченной совокупности слоев i (t), i=1,2,…,n (t) будем называть структурной моделью системы (t) и обозначать как g(t). В рамках сделанных предположений структурную модель g(t) назовем слоистой. Для простоты будем исходить из того, что слои i (t), i=1,2,…,n (t) не претерпевают структурных нарушений в области (t), а физические свойства материальной среды, как функции пространственных координат, в пределах слоя изменяются незначительно.
Пусть далее на временном отрезке [0,T] от момента зарождения системы (t) до сегодняшнего дня, может быть задана временная сетка:
t = {ti / ti = ti-1 +h ti, h ti 0, i=1,2,…,n t ;
t0 = 0, tnt = T}.
Узлы сетки t совпадают с моментами изменения внешних условий формирования системы, приводящими к структурным изменениям в ней, которые заключаются, например, либо во включении в нее нового слоя, либо в исключении из нее некоторого конечного количества уже имевшихся в системе слоев. На любом временном интервале Ti =(ti-1,ti ), i=1,2,…,n t в системе может происходить формирование только одного нового слоя. Свойства и пространственные границы всех слоев, входящих в систему, могут изменяться во времени.
взаимная диффузия различных по характеристик ам материальных сред предполагается пренебрежимо малой.
Зададим для определенности в области (t) декартовую1 систему координат OXYZ. Каждой пространственной точке A(x,y,z)(t), t[0,T] может быть поставлен в соответствие вектор (x,y,z,t)Es - вектор физических характеристик, определяющих локальные свойства системы (t) в этой точке в момент времени t. Здесь и далее Es – s-мерное евклидово пространство. s - количество указанных характеристик, часть из которых определяет свойства структурных элементов материальной среды, а часть свойства протекающих в системе процессов.
Будем считать справедливыми следующие утверждения:
1. Скорость формирования системы (t) и скорости протекающих в ней процессов очень малы по сравнению с процессами, проходящими в масштабе реального времени. Как следствие, временной промежуток [0,T] слишком велик для того, чтобы имелась возможность непосредственно наблюдать весь жизненный цикл системы от ее зарождения до настоящего времени. Можно считать, что система (t) потенциально доступна для непосредственного комплексного изучения (посредством прямых измерений) лишь в настоящий момент времени2 (при t=T).
2. Значения компонент вектора (x,y,z,T), как правило, недоступны для прямых измерений во всей области (T). Такие измерения возможны лишь для отдельных характеристик и лишь, вообще говоря, для конечного множества точек A(x,y,z)(T).
3. Существуют комплексные косвенные методы исследования, позволяющие с определенной точностью определить временной отрезок [0,T], задать временную сетку t, и восстановить на сетке t внешние условия и историю формирования системы (t). Как следствие t[0,T] имеется возможность построить структурную модель системы g(t) и определить для временных интервалов Ti =(ti-1,ti ), i=1,2,…,n t возможные свойства слоя i (t), i=1,2,…,n (t) при вхождении его в систему.
4. Процессы, протекающие в исследуемой системе, в достаточной мере изучены, что позволяет при необходимости использовать ряд упрощенных эмпирических законов, количественно описывающих изменение ее свойств с течением времени.
Для определенности будем считать, не умаляя общности, что пространственная область (t) здесь и далее задается в каждый момент времени t[0,T] условиями3 :
xmin x xmax, ymin x ymax, (2) 0 z zmax (x, y, t).
Перейдем к формулировке основной задачи. Предположим, что g(t) – структурная модель исследуемой системы известна для любого t[0,T]. Пусть имеется конечное множество пространственно локализованных подобластей Vk (T), k=0,1,2,…,n V, каждая из которых пересекается с одним и тем же множеством слоев структурной модели g(T). Для определенности будем считать, что подобласти Vk представляют собой вертикальные одномерные разрезы области (T).
Vk = { (xk,yk,,z)(T) / 0 z z max(Vk,T)}, k=0,1,2,…,n V.
Для каждого вертикального разреза Vk определена одномерная сетка:
kz = {zk,i / zk,i = zk,i-1 +h zk,i, h zk,i 0, i=1,2,…,n zk ;
zk,0 = 0, zk,nzk zmax(Vk,T)}.
Тип используемой системы координат зависит от специфик и пространственного строения исследуемой системы.
Впрочем, именно этот момент времени и имеет обычно наибольшее практическое значение.
Ось OZ считается направленной вниз.
Vk (T), Вертикальные разрезы k=1,2,…,n V, будем называть экспериментальными. Для них в узлах zk,i соответствующей сетки kz известны (в результате прямых или косвенных измерений) значения P* (xk,yk,zk,i,T) – значения целевой характеристики в настоящее время. Эти значения будем называть полевыми данными, а совокупность значений P* (xk,yk,zk,i,T), i=1,2,…,n zk обозначать как P*k (z) и называть полевым дискретным распределением целевой характеристики для разреза Vk.
Вертикальный разрез V0 (T) будем называть прогнозным. Для этого разреза полевые данные не заданы.
Основная задача, может быть теперь сформулирована следующим образом.
Задача P Пусть для каждого вертикального разреза Vk (T), k=1,2,…,n V известно полевое дискретное распределение P* k (z) целевой характеристики, заданное в узлах соответствующей сеткиkz.
Для заданного вертикального разреза V0 (T) требуется определить дискретное распределение P0 (z) в узлах сетки 0,z = {z0,i / z0,i = z0,i-1 +h z0,i, h z0,i 0, i=1,2,…,n z,0 ;
z0,0 = 0, z0,nz0 zmax(V0,T)}.
При разработке методологии решения Задачи P будем исходить из следующих предположений.
1. Получение требуемого результата посредством интерполяции известных полевых распределений P*k (z), k=1,2,…,n V целевой характеристики в прогнозный вертикальный разрез V0 не представляется возможным. Решение задачи должно происходить в контексте всего периода эволюционного развития системы.
2. Отсутствует достаточно полная и точная априорная количественная информация о физических свойствах структурных элементов системы в настоящее время. Количество n V экспериментальных вертикальных разрезов Vk относительно невелико.
3. Может быть построена (например, на основе тех или иных фундаментальных законов) физически интерпретируемая математическая модель, позволяющая находить P(x,y,z,T) - модельное распределение целевой характеристики в области (T) в контексте эволюционного развития системы. Как правило, практическое построение P(x,y,z,T) осуществляется в результате аппроксимации модели дискретной моделью (например, разностная аппроксимация краевой задачи для уравнения математической физики) с ее последующим численным исследованием на ЭВМ.
Представим итоговую зависимость P(x,y,z,T) от модельных параметров в явном виде:
P = G(), (3) где: P=P(T)PEn – модельное дискретное распределение целевой характеристики;
XEm - вектор варьируемых модельных параметров, соответствующих отдельным существенным для моделируемого процесса физическим характеристикам системы;
G() – оператор, осуществляющий отображение из X - пространства модельных параметров в P - пространство модельных распределений.
Задачу (3) будем называть задачей прямого моделирования или прямой задачей, оператор G() – оператором прямой задачи.
4. Непосредственное получение P0 (z) искомого дискретного распределения целевой характеристики для заданного вертикального разреза V0 (T) с помощью операторного соотношения (3) практически неосуществимо по нескольким причинам.
Во-первых, даже если предположить, что построение точной модели изучаемых процессов является выполнимой задачей, эта модель, как правило, оказывается для реальных природных систем чрезвычайно громоздкой и сложной не только для аналитического, но и для численного исследования.
Во-вторых, из-за объективного недостатка данных о фактических свойствах структурных элементов системы и протекающих в ней процессов, что исключает возможность получения приемлемых результатов численного моделирования изучаемых процессов даже при наличии точной модели.
5. В пределах каждого слоя i (t), i=1,2,…,n (t) структурной модели g(t) для физических свойств среды характерна непрерывная и, как правило, незначительная изменчивость в пространственном отношении.
В основу разрабатываемой в данной работе методологии решения Задачи P положены следующие принципы.
Принцип эволюционного развития Постановка и решение прямой задачи (3) должны осуществляться в контексте эволюционного развития системы (t) с учетом динамики изменения ее структуры и свойств образующей ее материальной среды.
Принцип эффективного моделирования Для моделирования исследуемого процесса должны использоваться эффективные модели. Под эффективной моделью понимается такая физически обоснованная модель исследуемого процесса, построение которой осуществляется для обеспечения адекватного моделирования заданной (целевой) характеристики. При этом перед эффективной моделью не ставится задача точного описания моделируемого процесса в целом. Адекватность моделирования понимается в том смысле, что эффективная модель должна быть по тенциально настраиваемой (калибруемой) на имеющиеся полевые данные о значениях целевой характеристики.
Другими словами, каковы бы ни были полевые данные, всегда можно выбрать такие значения модельных параметров (из множества их возможных значений), при которых моделируемые значения целевой характеристики будут в достаточной степени близки к полевым данным. Эффективная модель становится функционально адекватной лишь после ее калибровки (параметрической идентификации модели).
Принцип эффективного моделирования носит общий характер и имеет отношение не только математическим моделям, но и к любым другим моделям. По сравнению с точными моделями, эффективные модели, как правило, существенно более просты с точки зрения возможности их теоретического и численного анализа, что является их основным преимуществом. С другой стороны они имеют целенаправленный характер и не могут служить для полноценного изучения моделируемых процессов в целом.
Обязательным условием для рассматриваемых в работе эффективных моделей является возможность их естественной физической интерпретации.
Необходимость осуществления параметрической идентификации (калибровки) эффективной модели предполагает постановку и решение обратной задачи, которая в общем случае может быть сформулирована в следующем виде:
Задача Пусть в области (T) известно полевое дискретное распределение P* (x,y,z,T) целевой характеристики. Требуется найти такой вектор * X, для ко торого выполнено условие:
G( * ) = P* (4) Обратная задача (4), как правило, не является корректной и для ее решения требуется применение регуляризирующих подходов. В работе используется один из таких по дходов, ориентированный на построение так называемых -квазирешений.
Необходимо также обеспечить приемлемую вычислительную трудоемкость решения задачи (4). В этой связи используется следующий принцип.
Принцип минимальной сложности В соответствии с этим принципом используемые для решения поставленной задачи эффективные модели должны быть предельно просты с вычислительной точки зрения. В связи с этим из принятого класса эффективных моделей необходимо выбрать модель по возможности минимальной вычислительной сложности.
Реализация такой модельной идентификации в работе осуществляется посредством рациональной минимизации количества варьируемых модельных параметров прямой задачи за счет исключения из их числа тех параметров, влияние изменения значений которых на результат моделирования незначительно. С этой же целью осуществляется построение структурной модели минимально допустимой сложности за счет рационального уменьшения общего числа входящих в нее структурных элементов – слоев i (t), i=1,2,…,n (t). В результате достигается также большая сбалансированность сложности принятых модельных описаний и полевых данных.
Принцип пространственной локализации Физические свойства системы, ее структура могут объективно способствовать специфической пространственной ориентации исследуемого процесса. В этом случае, и сбор данных, и моделирование желательно осуществлять в соответствии с имеющей место пространственной спецификой. Кроме того, по объективным причинам полевые данные могут быть пространственно локализованы или пространственно ориентированы, что также должно по возможности учитываться при разработке вычислительных схем анализа этих данных.
В данной работе при формулировке Задачи P было сделано предположение, что полевые данные представлены только в вертикальных разрезах Vk (T), k=0,1,2,…,n V. Как следствие, соответствующие полевые дискретные распределения P*k (z), k=1,2,…,n V пространственно ориентированы вдоль оси OZ.
В этой связи естественно произвести декомпозицию исходной математической модели на одномерную математическую модель, ориентированную на описание «проекции» исследуемого процесса на ось OZ и двумерную – ориентированную на описание «проекции» процесса на горизонтальную плоскость OXY. Одномерную модель, полученную в результате такой декомпозиции, будем называть «вертикальной» моделью, а двумерную модель – «латеральной». При этом учет в вертикальной модели латеральных составляющих процесса осуществляется за счет введения в модель дополнительных модельных параметров. Аналогичным образом осуществляется учет вертикальных составляющих процесса в латеральной модели.
Предположим, что указанная декомпозиция проведена и построена вертикальная математическая модель исследуемого процесса. Аналогично прямой задаче (3) запишем вертикальную задачу в виде:
Pk(z)= Gz( (k)), k=1,2,…,n V (5) где: Pk(z)=Pk(z,T) – модельное распределение характеристики в вертикальном разрезе Vk (T), k=1,2,…,n V, а (k) - соответствующий вектор модельных параметров.
Задачу (5) будем называть одномерной задачей прямого моделирования или вертикальной прямой задачей, а оператор Gz() – оператором вертикальной прямой задачи.
Калибровку вертикальной модели естественно проводить отдельно для каждого вертикального разреза Vk (T), k=1,2,…,n V, решая обратную задачу в приведенной ниже постановке:
Задача z Пусть для вертикального разреза Vk (T) известно полевое дискретное распределение P*k (z), целевой характеристики.
Требуется найти такой вектор (k)* X, для которого выполнено условие:
Gz( (k)* ) = P*k (z) (6) В отношении задачи (6) могут быть приведены те же комментарии, что и к задаче (4) с точки зрения ее некорректности. Однако в вычислительном отношении нахождение численного решения задачи (6), безусловно, будет существенно менее трудоемким.
В контексте задачи (6) вертикальные разрезы Vk (T), k=1,2,…,n V будем называть также калибровочными разрезами, а вектора (k)* - соответствующими калибровочными векторами модельных параметров.
Общая схема методологии решения задачи численного прогнозирования характеристик природных динамических систем Представим общую схему (рис.1) предлагаемой методологии решения исследуемой задачи в виде конечного числа выполняемых в определенной последовательности этапов1 (все требуемые при выполнении этих э тапов действия предполагаются выполнимыми) и приведем необходимые комментарии к основным элементам (этапам) этой схемы.
1. Задание исходной информации о системе (t) и ее характеристиках. На этом этапе доступная априорная фактическая или экспертная информация о системе анализируется и подвергается необходимой предобработке с тем, чтобы появилась возможность задать все необходимые для исследования сведения. В том числе:
- временной отрезок [0,T];
- временную сетку t = {ti / ti = ti-1 +h ti, h ti 0, i=1,2,…,n t ;
t0 =0, tnt =T};
- внешние условия формирования системы для каждого из временных интервалов Ti =(ti-1,ti ), i=1,2,…,n t ;
- количество слоев, определяющих структурную модель системы;
- возможные свойства слоя i (t), i=1,2,…,n (t) при вхождении его в исследуемую систему;
- конечное множество калибровочных разрезов Vk и P*k (z) - имеющих место в настоящее время (t=T) дискретных полевых распределений значений целевой характеристики для каждого Vk, k=1,2,…,n V.
2. Построение g(T) - структурной модели системы осуществляется экспертно на основании информации, заданной на первом этапе.
3. Формирование эффективной математической модели исследуемого процесса и осуществление при необходимости ее декомпозиции. Определение перечня и задание диапазона возможных изменений значений модельных параметров. Этот этап полностью зависит от специфических особенностей конкретной системы и исследуемого процесса, а также структуры полевых данных. Построение модели осуществляется в масштабе времени развития системы.
первый и второй этапы общей схемы, будучи обязательными эл ементами методологии, не являются предметом исследов ания в данной работе (исходные данные и структурную модель системы будем считать заданными). То же относится и к этапу 9.
Сбор и предобработка 1 Модельная исходной информации о идентификация.
системе (t) и ее Снижение количества характеристиках модельных параметров Построение g(t) - 2 Параметрическая структурной модели идентификация.
системы Калибровка модели по полевым данным Построение 3 Построение эффективной модели прогноза в заданной процесса. пространственной области Прямая задача. 4 Анализ Численные методы полученных результатов решения Обратная задача. 5 Уточнение прогноза Численные методы при поступлении решения дополнительной информации - блоки, выполняемые при разработке вычислительной технологии;
- блоки, выполняемые в интерактивном режиме, с участием человека;
- блоки, выполняемые в автоматическом режиме;
- не обязательно выполняемые блоки;
- завершение работы.
Рис. 1. Общая схема методологии численного прогнозирования характеристик природных динамических систем 4. Постановка прямой задачи осуществляется в соответствии с выбранной эффективной математической моделью процесса. Построение численных методов решения прямой задачи должно учитывать ее математическую специфику.
5. Постановка обратной задачи не зависит от специфики исследуемой системы и может быть осуществлена в общем случае. Для решения обратной задачи предлагается использовать, как отмечалось выше, регуляризирующий подход, ориентированный на построение так называемых -квазирешений.
Рассмотренные первые пять этапов общей схемы являются подготовительными этапами, предваряющими собственно процесс решения Задачи Р. Все они выполняются при непосредственном участии специалистов в соответствующей предметной области. При этом третий, четвертый и пятый этапы выполняются, как правило, лишь один раз при подготовке к численному решению задачи. При практическом применении методологии эти этапы могут быть исключены из общей схемы.
6. Упрощение используемых эффективных моделей (модельная идентификация) осуществляется посредством уменьшения количества модельных параметров в результате исключения из числа варьируемых модельных параметров тех, изменение значений которых не оказывает существенного влияния на качество калибровки модели исследуемого процесса. Аналогичным образом исследуется возможность объединения (в один слой) любых дву х соседних слоев i (t), i+1 (t), i=1,2,…,n (t)-1 в структурной модели системы g (t). Весь процесс организовывается, как правило, в интерактивном режиме. В том числе и задание необходимых критериев исключения модельных параметров. После завершения этого этапа формируется окончательный (рабочий) вариант структурной модели системы и определяется рациональное количество варьируемых модельных параметров прямой задачи.
7. Калибровка (параметрическая идентификация) вертикальной математической модели для каждого калибровочного разреза Vk, k=1,2,…,n V осуществляется посредством решения обратной задачи (6). В результате калибровки определяются значения компонент калибровочных векторов (k)*. В зависимости от конкретных целей, калибровка может проводиться, либо в автономном режиме для каждого калибровочного разреза, либо в режиме одновременной калибровки всех калибровочных разрезов с дополнительными условиями на согласование получаемых результатов (региональная калибровка). Установка значений управляющих процессом констант осуществляется в интерактивном режиме.
8. Построение прогноза дискретного распределения P0 (z) в заданной пространственной области V0. Этот этап выполняется полностью в автоматическом режиме и состоит из дву х шагов.
На первом шаге посредством послойной интерполяции значений компонент калибровочных векторов (k)*, k=1,2,…,n V в вертикальный разрез V0 определяются значения компонент вектора (0)* - вектора модельных параметров для разреза V0.
На втором шаге в результате решения прямой задачи:
P0 (z)= Gz( (0)* ) осуществляется построение требуемого прогноза.
9. Анализ полученных результатов прогнозирования осуществляется экспертно в интерактивном режиме. При неудовлетворительном заключении может производиться дополнительный анализ данных о системе и ее свойствах и целесообразная корректировка структурной модели g (t) с целью повторения процесса прогнозирования, начиная с этапа 7.
10. Этап уточнения прогноза является важным дополнительным элементом предлагаемой методологии. Его выполнение актуально при появлении различного рода дополнительной информации о свойствах исследуемой системы в окрестности разреза V0 в случае значимого расхождения прогнозных и наблюдаемых значений целевой характеристики. Один из возможных подходов к построению процедуры такого уточнения предложен в главе 5 данной работы.
Рассмотренная общая схема методологии численного прогнозирования характеристик природной динамической системы имеет общий характер. Она позволяет в полной мере реализовать сформулированную концепцию общего подхода к решению Задачи P. Вместе с тем, очевидно, что ее практическая реализация возможна лишь применительно к конкретным природным динамическим системам.
В последующих главах работы на основе этой общей схемы осуществляется разработка методологии численного прогнозирования ГФД в осадочных бассейнах земли. Важно отметить, что в контексте Задачи Р все предлагаемые при этом модели и алгоритмы носят достаточно общий характер и могут быть использованы при исследовании любых природных систем, удовлетворяющих сформулированным выше предположениям об их свойствах, за исключением учитывающей специфику ГФС постановки прямой задачи и предлагаемых методов ее решения (Главы 2 и 3).
Во второй главе осуществляется постановка задачи прогнозирования ГФД в осадочных бассейнах земли и формулируются эффективные модели ГФС с целью получения прогноза эксцесса ГФД (превышения ГФД над гидростатическим давлением) в заданной подобласти осадочного бассейна. В частности, в области бурения новой скважины.
Наличие такого прогноза помогает избежать аварий на буровых установках, снизить риск потери больших материальных средств и дорогостоящего оборудования.
Вместе с тем практически используемые в настоящее время классические подходы к прогнозированию ГФД основаны, как правило, на упрощенных эмпирических законах, ч то нередко приводит к серьезным ошибкам при сложной истории развития осадочного бассейна, не учитываемой в упрощенных моделях.
В контексте рассмотренных ранее общих построений и предположений и при использовании прежних обозначений задача прогнозирования ГФД может быть сформулирована аналогично Задаче Р.
Основная задача Пусть в некоторой пространственной области1 (T), принадлежащей осадочному бассейну, задано конечное множество Vk (T), k=0,1,2,…,n V, - скважин (вертикальных одномерных разрезов области (T)). Для каждой скважины Vk в узлах область (T) по-прежнему определяется условиями (2) в выбранной декартовой системе координат соответствующей одномерной сетки kz известны значения P* (xk,yk,zk,i,T), i=1,2,…,n zk – значения ГФД (полевые данные), образующие в совокупности вектор P*k (z), задающий полевое дискретное распределение ГФД для скважины Vk.
Требуется определить P*0 (z) – прогноз дискретного распределения ГФД на сетке 0z для прогнозной скважины (вертикального разреза) V0 (T).
Существенно, что рассматривается проблема формирования прогноза для целей разведочного бурения. Это позволяет считать, что естественный эволюционный характер развития осадочного бассейна полностью определяет текущее состояние системы, поскольку воздействие человеческой деятельности на него на этом этапе исследования бассейна пренебрежимо мало.
Построение необходимых для исследования моделей должно осуществляться в геологическом масштабе времени, что, как отмечалось ранее, позволяет отнести рассматриваемую задачу к задачам бассейнового моделирования. Работы в этой области активно ведутся уже несколько десятков лет. В результате в настоящее время имеются развитые инструменты бассейнового моделирования, основанные, в том числе, на численных решениях многофазных задач флюидодинамики. Вопросы моделирования и прогнозирования ГФД с учетом истории геологического развития ГФС исследовались, в частности, в работах Буряковского Л.А., Джафарова И.С., Джеваншира Р.Д. и других авторов. В работах Lеrchе I., Zhao K., Yu Z и других разрабатывались подходы к прогнозированию характеристик ГФС на основе постановки и решения обратных задач. Однако, построение реализуемой в реальном масштабе времени и практически эффективной методологии численного прогнозирования ГФД остается актуальным.
Основным объектом моделирования при разработке методологии численного прогнозировании ГФД является ГФС. В соответствии с общей вычислительной схемой методологии численного прогнозирования характеристик природных систем необходимо, прежде всего, сформировать структурную модель или в данном случае геологическую модель ГФС и математическую модель исследуемого ФД процесса.
В содержательном плане построение геологической модели ГФС сводится к решению проблемы нелокального выделения геологических объектов, которая исследовалась в работах Губермана Ш.А., Карогодина Ю.Н., Мушина И.А. и целого ряда других авторов.
В данной работе при построении геологической модели ГФС в качестве достаточного уровня детальности для целей проводимого исследования принят уровень формационных ассоциаций (формаций) геологических объектов. Способы и методы выделения формационных элементов геологической модели не являются предметом обсуждения в данной работе. Отметим лишь, что для этих целей могут быть использованы известные в геологии и геофизике методы структурно-формационной интерпретации (СФИ) комплекса геолого-геофизической информации.
Применительно к формациям, как структурным элементам геологической модели осадочного бассейна, оправданно предположение о пространственно-временной упорядоченности и квазиоднородности формаций (слоев). Наиболее древние слои находятся в осадочном бассейне, как правило, глубже менее древних. Физические свойства осадочных пород в пределах формации мало изменчивы в пространственном отношении. Поэтому будем предполагать, что в вертикальном направлении в любой момент времени физические свойства осадочных пород, образующих формацию, постоянны.
Рис.2. Слоистая структура геологической модели осадочного бассейна (фрагмент разреза) Начальная геологическая модель (рис.2) осадочного бассейна может содержать в общем случае до сотни различных формаций осадочных пород. С позиций эффективного подхо да к моделированию предполагается, что каждая из выделенных формаций может быть отнесена по характерным свойствам (составу, структуре и т.д.) образующих их осадочных пород к одному из трех базовых типов (глинистому, песчанистому и карбонатному).
Существенно, что для большинства осадочных бассейнов по оценкам специалистов лишь 5-10% исследуемого интервала по глубине бассейна может обеспечивать значимый уровень латерального (в плоскости ХОУ) оттока флюида. Эти 5-10% глубинного интервала представлены формациями, которые будем называть латерально-проводящими каналами. Таким образом, можно считать, что почти на всем глубинном интервале осадочного бассейна фильтрация флюида происходит в вертикальном направлении.
При построении геологической модели учитываются также различного рода структурные нарушения. В том числе вызванные возможными перерывами в осадконакоплении - имевшими место в истории формирования осадочного бассейна периодами подъема части уже сформировавшихся формаций выше уровня моря и их полным или частичным выветриванием (эрозией). К структурным нарушениям относятся также наблюдаемые в осадочном бассейне в настоящее время разломы (рис.3) – тектонические разрывы и смещения в пространственном расположении формаций.
Предполагается, что на основе содержательной обработки и интерпретации всего комплекса доступной геолого-геофизической информации об осадочном бассейне может быть определена с некоторой степенью погрешности и другая необходимая для проведения дальнейших исследований информация. В том числе возраст и период начального формирования, предполагаемый темп осадконакопления и диапазоны возможных значений физических характеристик среды для каждой учитываемой в геологической модели формации, включая подвергшиеся полной или частичной эрозии. Другими словами может быть восстановлена в эффективном смысле вся история формирования исследуемого осадочного бассейна с учетом возможных перерывов в осадконакоплении. В результате для временного отрезка [0,T] может быть задана временная сетка:
t = {ti / ti = ti-1 +h ti, h ti 0, i=1,2,…,n t ;
t0 = 0, tnt = T}.
Узлы сетки t совпадают с моментами структурных изменений в системе. Как и в общем случае, эти изменения заключаются, например, либо в начале процесса вхождения в систему нового слоя (формации), либо в исключении из нее некоторого конечного количества уже имевшихся в системе слоев, а также в возникновении разломов, либо в некоторой комбинации указанных событий. На любом временном интервале Ti =(ti-1,ti ), i=1,2,…,n t может происходить процесс формирования в системе км ] Гл уб ин а [ Расстояние [10 * км] b a Рис. 3. а – разрез фрагмента геологической модели с примерами разломов формаций b – вид сверху по верхней границе третьей снизу формации. Линия разреза нанесена пунктиром.
только одного нового слоя.
Процесс геологического формирования осадочного бассейна - от его зарождения до настоящего времени – сопровождается происходящими в нем процессами различного характера, в результате которых свойства и пространственные границы формаций, образующих систему изменяются. Одним из важнейших среди них является ФД процесс, который заключается в миграции флюида (жидкости и газа) в толщах осадочных пород под воздействием поля ГФД.
Теоретические основы процессов уплотнения осадочных пород, перехода углеводородов при нагревании материнских пород (осадочных толщ, содержащих углеводороды) из твердой в жидкую и газообразную фазы (УВ генерация), миграции флюидов в осадочных толщах и т.д. рассматривали в своих работах Александров Б.Л., Гуревич А.Е., Добрынин В.А., Жузе Т.П., Коновалов А.Н., Серебряков В.А., Шестаков В.М., Bachmat Y., Bear J., Magara K., Tissot B.P., Verwei J.M., Welte D.H. и другие российские и зарубежные авторы. Все эти процессы имеют большое научное и практическое значение с точки зрения возникновения в ГФС зон сверхгидростатических ГФД и в эффективном смысле учитываются при разработке методологии численного прогнозирования ГФД.
В данной работе для построения эффективной математической модели реального многофазного ФД процесса предлагается использовать уравнение фильтрации однофазного флюида - несжимаемой жидкости. Такое существенное упрощение носит принципиальный характер и объясняется необходимостью построения как можно более простой эффективной модели с целью обеспечения возможности решения Основной задачи в приемлемое время. Все рассмотрения проводятся в принадлежащей осадочному бассейну трехмерной области (t), определяемой условиями (2), на временном промежутке [0,T]. Дневная поверхность области (t) (при z=0) совпадает с уровнем моря. Построение математической модели осуществляется на основе фундаментального закона сохранения (консервации) масс.
Скорость фильтрации флюида определяется линейным законом Дарси. В общем виде предлагаемая математическая модель представляет собой линейное параболическое уравнение математической физики с переменными коэффициентами:
A(P/t) = (k x(P/x))/x + (k y (P/y))/y+ (k z(P/z))/z + B, (7) где: A(x,y,z,t), B(x,y,z,t) – формально функции пространственных координат и времени, а по существу - свойств осадочной породы и флюида. k x(x,y,z,t), ky (x,y,z,t), k z(x,y,z,t) – коэффициенты проводимости среды, определяющие ее фильтрационные свойства.
P(x,y,z,t) – превышение (эксцесс) ГФД над гидростатическим давлением – моделируемая характеристика ГФС.
Сформулирована краевая задача для уравнения (7). Однако с практической точки зрения ее решение для целей численного моделирования пространственного распределения ГФД в области (T) не представляется целесообразным ввиду недостаточной обеспеченности априорными данными с одной стороны и большой вычислительной трудоемкости с другой.
В третьей главе рассматриваются вопросы декомпозиции предложенной трехмерной (по пространственным переменным) эффективной модели флюидодинамического процесса на одномерную и двухмерную модели.
Формулируются соответствующие прямые задачи, осуществляется построение методов решения сформулированных прямых задач.
С учетом указанных ранее особенностей протекания ФД процесса в осадочном бассейне, вместо модели (7) вводятся в рассмотрение одномерная (вертикальная) и двумерная (латеральная) модели:
A(P/t) = (k z(P/z))/z + q L + B, (8) z A(P/t) = (k x(P/x))/x+ (k y (P/y))/y + q + B, (9) L Здесь модельный параметр q предназначен для эффективной компенсации не учитываемой в уравнении (8) латеральной составляющей потока флюида, а модельный параметр q z – для эффективной компенсации не учитываемой в уравнении (9) вертикальной составляющей потока.
В результате, «вертикальная» 1 прямая задача формулируется в виде краевой задачи для уравнения (8):
A(P/t) = (k z(P/z))/z + q L + B, t[0,T], z[0, z max(x, y, t)], P(z)=0 z=0, (10) (P/z) z=zmax = 0, P(z) = P0 (z) t=0.
Необходимо подчеркнуть, что zmax(x, y, t) - нижняя граница пространственной области решения задачи (10) с течением времени изменяется. Как правило, она увеличивается – разрез становится глубже. Однако в ряде случаев может и уменьшаться на некоторый период времени. Например, при относительном снижении уровня моря.
Пусть L (t)(t) – подобласть области (t), соответствующая L-му латерально-проводящему каналу, и GL = GL0 GL1 GL2 – боковая граница L (t).
Соответствующая «латеральная» прямая задача формулируется как краевая задача для уравнения (9). В несколько в упрощенном виде она приведена ниже:
A(P/t) = (k x(P/x))/x+ (k y (P/y))/y + q z + B, t[ tL,T], x[ x min, xmax], y[ ymin, ymax], L P/n = 0 (x,y,zL )G 0, (11) L P/n = P/n (x,y,zL )G 1, L P = const (x,y,zL )G 2, P = P0 t=tL.
Здесь: tL – время начала процесса образования латерально-проводящего канала;
zL – глубина залегания канала в осадочном бассейне;
GL0 – часть границы области L (t), поток флюида через которую отсутствует;
GL1 – часть границы, через которую происходит поток флюида в «область разгрузки»;
n - усредненное расстояние от границы GL1 области L (t) до области разгрузки2 ;
P – усредненное значение разности эксцессов ГФД на GL1 и в области разгрузки;
GL2 – часть границы, на которой эксцесс ГФД со временем не изменяется;
n – нормальный вектор к соответствующей части боковой границы области L (t);
P0 – начальное пространственное распределение эксцесса ГФД в области L (tL ).
В общем случае zL зависит не только от времени, но и (см., например, рис.3) от пространственных координат x и y. Кроме того, мощность (расстояние от нижней до верхней границы) латерально-проводящего канала также зависит от x, y и t. В приведенной формулировке задачи (11) это в явном виде не отражено, однако, при разработке метода ее решения указанные особенности задачи учитываются алгоритмически.
Вопросы решения краевых задач математической физики, задач фильтрации однофазных и многофазных флюидов в сплошных однородных и неоднородных средах рассматривали в своих работах Бахвалов Н.С., Бер Я., Бобков В.В., Годунов значения пространственных координат x и y фиксированы Под областью разгрузки понимается некоторая область осадочного бассейна с пониженным значением эксцесса ГФД, с которой у L имеется флюидодинамическая связь.
С.К., Коновалов А.Н., Лейбензон Л.С., Марчук Г.И., Самарский А.А., Флетчер К., Шестаков В.М. и многие другие российские и зарубежные авторы.
В данной работе решение вертикальной и латеральной прямых задач (10) и (11) осуществляется на основе классических сеточных методов. Основная особенность предлагаемых подходов состоит в том, что процесс решения на отрезке [0,T] разбивается на конечное число последовательно выполняемых макрошагов. Каждый макрошаг соответствует временному интервалу Ti =(ti-1,ti ), i=1,2,…,n t временной сетки t (продолжительность i-го макрошага равна Ti = ti - ti-1 ). Для вертикальной задачи, таким образом, число макрошагов (см. рис.4) равно n t, а для латеральной задачи количеству временных интервалов Ti на отрезке [tL,T] с момента образования латерального канала.
На каждом макрошаге решение осуществляется в два этапа.
Для конкретизации рассмотрим макрошаг, соответствующий временному интервалу Ti, на котором в систему добавляется новая формация.
Первый этап (э тап нагрузки).
t0 t1 t2 t3 t4 t5 t6 T t Предполагается, что все изменения в геологической модели, обусловленные формированием нового слоя, процессами погружения и уплотнения пород, их нагреванием, возможной генерацией углеводородов и т.п., произошедшие за временной промежуток Ti, полностью завершены.
ГФС находится в “нагруженном” состоянии. Исхо дя из этого предположения, корректируются глубина залегания, мощности и свойства всех формаций, образующих к 1 2 3 45 6 номера макрошагов моменту времени t=ti геологическую z модель пространственной области Рис. 4 Схема разбиения процесса решения решения. Корректировка производится прямой задачи на макрошаги для на основе тех или иных эмпирических вертикального разреза из семи формаций.
методов в зависимости от типа конкретных формаций. Затем для каждой формации, входящей в откорректированную геологическую модель, производится расчет текущих значений коэффициентов уравнений (8) или (9). На основании пространственного распределения эксцесса ГФД, полученного на предыдущем макрошаге (при t=ti-1 ), и дополнительной нагрузки, которую получила система к моменту времени t=ti, формируется P0 - начальное пространственное распределение эксцесса ГФД для решения краевой задачи на временном интервале Ti.
Второй этап (этап разгрузки). Геологическая модель и свойства ее структурных элементов, определенные на первом этапе, считаются неизменными.
Посредством численного решения прямой задачи (10) или (11) на временном интервале Ti в текущей пространственной области осуществляется «разгрузка» системы за счет ФД процесса. После окончания решения прямой задачи получаем пространственное распределение эксцесса ГФД в области решения на момент завершения макрошага.
Для решения прямой задачи (10) на втором этапе каждого макрошага рассмотренной вычислительной схемы применяется дву хслойная неявная четырехточечная разностная схема. Для решения прямой задачи (11) – метод, в основе которого лежит одна из разностных схем метода переменных направлений, называемая продольно-поперечной разностной схемой Писмена-Рэчфорда. В работе обсуждаются имеющиеся особенности в реализации названных методов, обусловленные, например, имеющей место в общем случае (см. рис.3) пространственной изменчивостью глубины залегания и мощности латерального канала.
Рассмотренный подход к решению прямых задач в целом является эффективным. К нему не предъявляется высоких требований по точности получаемых решений. Это обусловлено невысокой, как правило, точностью исходных данных, а также эффективностью используемых моделей. Главное требование, предъявляемое к методам решения – вычислительная устойчивость, что полностью обеспечено в предлагаемой вычислительной схеме. Существенное влияние на результаты моделирования оказывает выбор тех или иных значений модельных параметров.
Задание множества модельных параметров является важной проблемой. К их числу необходимо отнести те свойства среды, которые с одной стороны являются существенными для моделируемого процесса, а с другой стороны, количественная оценка значений которых для исследуемой области затруднительна (известен, например, лишь диапазон их возможных значений). Вместе с тем, общее количество модельных параметров не должно быть слишком большим. В данной работе для каждого формационного элемента геологической модели выделяется пять модельных параметров. Кроме уже упомянутых ранее параметров qL (для вертикальной задачи) и q z (для латеральной задачи), это начальная пористость осадочной породы, коэффициент ее уплотнения, константа удельной поверхности пор, начальное значение УВ-потенциала (характеризует потенциальную способность осадочной породы к генерации нефти или газа). В общем случае в латеральной задаче модельный параметр qz вводится отдельно для каждой пересекающей данный латерально-проводящий канал (см. рис.6) скважины Vk (T), k=1,2,…,n V.
Пусть - вектор модельных параметров для вертикальной задачи.
=(11,21,…,r1,1 2,22,…,r2,…,1 Nf,2 Nf,…,r Nf)T (12) NF – общее количество формаций, выделенных в геологической модели области, r – количество модельных параметров для одной формации. Значения модельных параметров могут выбираться из интервалов возможных значений:
k i min ki ki max, k=1,2,…,r ;
i=1,2,…, NF. (13) X X, если компоненты вектора удовлетворяют условиям (13). Здесь X – конечномерное евклидово пространство векторов модельных параметров, X0 – ограниченное замкнутое подмножество X, определяемое условиями (13).
Пусть Gz() – оператор, осуществляющий, в результате решения вертикальной прямой задачи (10) с применением предлагаемой вычислительной схемы, однозначное отображение из X0 в P0 P, где P – конечномерное евклидово пространство векторов модельных данных, а P0 – ограниченное замкнутое подмножество P. В отношении оператора Gz() можно сказать, что он является результирующей эффективной численной моделью моделируемой характеристики.
Аналогично определяется оператор GL () для латеральной прямой задачи (11).
Тогда для вертикальной прямой задачи будем использовать запись в форме:
P = Gz(), X0 X, PP0 P. (14) Аналогично для латеральной задачи будет использоваться запись в форме:
PL = GL (L ), L XL0 XL, PL PL0 PL. (15) В ряде случаев может оказаться полезной Gz(), tTi предлагаемая в работе вычислительная схема совместного решения вертикальной и латеральной прямых задач. По существу это qL qz пример декомпозиционного подхода к решению трехмерной краевой задачи для уравнения (7).
GL (L ), tTi Схема реализуется в рамках уже рассмотренного выше общего подхода с разбиением временного отрезка [0,T] на Рис.5 Схема декомпозиционного конечное число макрошагов и разделения подхода.
флюидодинамического процесса на каждом шаге на этап нагрузки и э тап разгрузки. Отличие состоит в реализации на каждом шаге второго этапа - этапа разгрузки ГФС.
Основная мотивация организации совместного решения задач (14) и ( заключается в том, что значение эффективного модельного параметра qL для вертикальной прямой задачи может быть получено из решения соответствующей латеральной задачи, а значение эффективного модельного z параметра q для латеральной прямой задачи – из решения вертикальной задачи.
В результате этап разгрузки на каждом Рис.6. Схематичное изображение макрошаге выполняется в виде пространственной области, в которой итерационной процедуры (рис.5) осуществляется совместное решение поочередного решения вертикальных и вертикальных и латеральных задач.
латеральных прямых задач для всех Калибровочные скважины изображены в калибровочных скважин и латеральных виде вертикальных столбцов, каналов, представленных (рис.6) в области латеральные каналы – в виде трехмерных моделирования. Итерации завершаются, поверхностей.
размерность XL и PL в (15) отличаются в общем случае от размерности X и P в (14).
когда корректировки в значениях q L и q z становятся достаточно малыми.
В работе обсуждаются преимущества, общее построение и алгоритмические особенности совместного решения вертикальной и латеральной прямых задач.
В четвертой главе рассматривается уточненная постановка обратной задачи и методы ее решения, а также методы упрощения (модельной идентификации) используемых эффективных моделей. Все построения в главе носят общий характер и могут быть использованы при исследовании любых природных систем в рамках предложенной в Главе 1 общей методологии.
Приемлемое качество прямого моделирования будет обеспечено при условии достижения необходимой близости модельных и полевых данных в выбранной метрике. В условиях недостатка достоверной информации об истинных характеристиках моделируемого процесса и эффективности используемых моделей требуемый результат может быть достигнут при надлежащем выборе значений модельных параметров модели. Такой выбор будем называть параметрической идентификацией модели или калибровкой. Калибровка модели может быть осуществлена в результате постановки и решения соответствующей обратной задачи.
Пусть для калибровочной скважины Vk (T) известно P* – полевое дискретное распределение эксцесса ГФД. Тогда аналогично задаче (6) формальная постановка обратной задачи для вертикальной краевой задачи (10) может быть представлена в виде:
Найти вектор * X0, для которого выполнено условие:
Gz( * ) = P* (16) В общем случае задача (16) некорректна. Наиболее существенным обстоятельством с практической точки зрения является не единственность ее решения.
Обратные задачи часто возникают в естественнонаучных исследованиях и вопросы их практического решения чрезвычайно актуальны. Здесь следует, прежде всего, отметить основополагающие работы А.Н.Тихонова по теории регуляризирующих методов решения некорректных задач. Проблемы решения некорректных задач рассматривали в своих работах О.М.Алифанов, В.Я.Арсенин, А.Б.Бакушинский, А.В.Гончарский, В.К.Иванов, М.М.Лаврентьев, В.Н.Страхов и многие другие отечественные и зарубежные ученые.
В данной работе в качестве регуляризирующего подхода к решению обратной задачи (16) предлагается использовать, как уже отмечалось ранее, известную концепцию построения так называемых -квазирешений. Этот подхо д реализован в соответствии с общими положениями теории построения регуляризирующих алгоритмов решения обратных задач в рамках схемы метода наименьших квадратов.
Введем в рассмотрение векторную функцию:
R() = D( Gz() – P* )/ P*, (17) * где: D – диагональная матрица весовых коэффициентов;
P - евклидова норма вектора P*.
Определим функцию рассогласования модельных и полевых данных в виде:
F() = R() 2 (18) В качестве -квазирешения задачи (16) определим любой вектор, удовлетворяющий условию;
F( ), X0, (19) где: 0 – заданное число.
Нахождение -квазирешения задачи (16) предлагается осуществлять с помощью сходящейся итерационной процедуры вида:
(i+1) = (i)+qi S(i), S(i)X, (i)X0, i = 0,1,2,…, (20) (0) 0 (i) где: X – начальное приближение;
S –направление спуска для функции рассогласования F();
q i 0 длина шага в направлении S(i).
Проблемы нелинейной оптимизации, методы решения линейных и нелинейных систем являлись и являются предметом многочисленных исследований.
Н.С.Бахвалов, В.В.Воеводин, Дж.Голуб, Дж. Деммель, Дж.Дэннис, В.Н.Кублановская, А.А.Самарский,, Ч.Лоусон, Л.С.Лэсдон, Н.Н.Моисеев, Р.Шнабель – лишь неполный перечень авторов, работы которых оказали влияние на результаты данного исследования.
В данной работе предложены алгоритмы реализации процесса (20) для построения -квазирешения задачи (16) на основе методов градиентного поиска и методов типа метода Гаусса-Ньютона. В последнем случае итерационный процесс (20) реализуется для решения (как правило, переопределенной) нелинейной системы уравнений вида:
R() =0, 0P, X0, (21) где: 0 – нулевой вектор.
Предлагаемый метод решения (21) строит, начиная (0)X0, последовательность точек (0), (1),..., (i),..., принадлежащих X0, по правилу:
J( (i)) S(i) = - R( (i)),, i = 0,1,2,…, (22) (i+1) = (i)+qi S(i), S(i)X, (i)X0, q i E1.
где: J() - матрицa Якоби векторной функции R(), вычисленная в точке.
В работе исследуются вопросы обеспечения сходимости метода решения обратной задачи, возможности учета дополнительной априорной информации, особенности алгоритмической реализации различных элементов предлагаемого подхода. В том числе – проблемы определения длины шага qi, обеспечения выполнения условия (i)X0, вычисления компонент вектора градиента функции F() и элементов матрицы J(). С целью повышения вычислительной устойчивости определения направления S(i), в результате решения линейной системы в методе (22), предлагается осуществлять сингулярное разложение матрицы J() с последующим сингулярным анализом. Многочисленные вычислительные эксперименты на реальных и синтетических данных подтверждают вычислительную эффективность разработанных алгоритмов решения обратной задачи.
Предложенный метод позволяет осуществлять калибровку вертикальной прямой задачи на основе полевых данных для конкретной калибровочной скважины с целью получения соответствующего калибровочного вектора модельных параметров.
Метод носит общий характер и может быть использован также для калибровки латеральной задачи и любых других эффективных моделей, если соответствующая обратная задача представима в форме, аналогичной (16).
Без принципиальных изменений может быть поставлена и решена задача одновременной калибровки всех калибровочных скважин, имеющихся в исследуемой области осадочного бассейна. Однако в этом случае существенно возрастает размерность задачи, что отрицательно сказывается на практической эффективности метода. Возможные подходы к организации процедур согласованной региональной (многоскважинной) калибровки рассматриваются в Главе 5.
Существенной проблемой на предварительном этапе исследования является уменьшение общего количества варьируемых модельных параметров исследуемой модели, что способствует снижению вычислительной трудоемкости ее последующей калибровки. Поскольку изменение количества варьируемых модельных параметров фактически изменяет саму модель, речь идет об идентификации модели в выбранном классе моделей. Рассматривается два подхода к решению этой проблемы.
Один из них заключается в эффективной корректировке геологической модели с целью уменьшения количества входящих в нее формационных элементов, что осуществляется в результате проверки возможности объединения в одну формацию дву х соседних в пространственном отношении и относящихся к одному типу формаций. Исключение формации из геологической модели автоматически приводит к исключению из рассмотрения всех относящихся к этой формации модельных параметров. Критерием возможности такого объединения служит сохранение допустимого качества калибровки модели прямой задачи, которое оценивается, например, проверкой возможности выполнения условия:
F(), X0, (23) где: 0 назначаемый исследователем допустимый уровень рассогласования между модельными и полевыми данными.
Второй подхо д ориентирован на проверку возможности исключения из числа варьируемых модельных параметров того или иного параметра индивидуально, без корректировки геологической модели. Критерий возможности такого исключения прежний – сохранение допустимого качества калибровки.
Большое значение в процессе идентификации модели имеет анализ так называемой чувствительности модельных параметров. Количественная оценка чувствительности модельного параметра определяется, в частности, на основе оценки относительного изменения значения F() - функции рассогласования при некотором изменении значения модельного параметра. Чем больше эта оценка, тем более чувствительна модель к изменению данного модельного параметра - параметр более чувствителен. Вследствие нелинейности оператора прямого моделирования, эти оценки, в общем случае, имеют локальный характер. В работе формулируются основные требования и правила формирования оценок чувствительности параметров, соблюдение которых позволяет достаточно обоснованно и эффективно использовать их в процедурах калибровки и идентификации моделей.
Калибровка модели, оценка и анализ чувствительности ее параметров положены в основу предлагаемых алгоритмических процедур идентификации модели. Выполнение этих процедур в реальном вычислительном процессе осуществляется, как правило, в интерактивном режиме. В результате оказывается возможным упростить геологическую модель и уменьшить количество модельных параметров прямой задачи, что способствует повышению практической эффективности методологии в целом.
В пятой главе рассмотрены алгоритмические процедуры согласованной (региональной) калибровки модели в исследуемой области и формирования прогноза эксцесса ГФД для заданной прогнозной скважины, предложена процедура уточнения прогноза в процессе бурения прогнозной скважины. Все построения в пятой главе носят достаточно общий характер и могут быть использованы при исследовании и других природных систем в рамках общей схемы методологии численного прогнозирования характеристик природных динамических систем (Глава 1).
Разработка завершающей процедуры предлагаемой методологии численного прогнозирования ГФД – процедуры формирования прогноза распределения эксцесса ГФД базируется, в частности, на предположении, что при отсутствии в исследуемой области осадочного бассейна структурных нарушений, калибровочные значения модельных параметров эффективной численной модели процесса должны иметь в ней достаточно гладкое пространственное распределение с незначительной изменчивостью. Исходя из э того предположения, калибровочные значения модельных параметров могут быть проинтерполированы в область прогнозирования ГФД. Используя введенные ранее обозначения, представим схему формирования прогноза распределения эксцесса ГФД для прогнозной скважины в виде следующей последовательности действий.
1. Определить калибровочные векторы (k)* X0 для всех калибровочных скважин Vk (T), k=0,1,2,…,n V.
2. Посредством интерполяции калибровочных значений одноименных модельных параметров определить компоненты (0)X0 - вектора модельных данных для прогнозной скважины V0 (T).
3. Определить P* 0 (z) - дискретное распределение эксцесса ГФД (прогноз) для скважины V0 как решение вертикальной прямой задачи:
P*0 (z) = Gz ( (0)), (0)X0 X, PP0 P. (24) Приведенная схема нуждается в некоторых комментариях.
Множество калибровочных векторов (k)* X0, k=0,1,2,…,n V должно быть согласовано между собой. Вследствие не единственности решения обратной задачи гарантировать это согласование при индивидуальной калибровке каждой отдельной калибровочной скважины невозможно без введения дополнительных априорных условий. В качестве основного априорного условия предлагается использовать условие минимальной вариабельности значений одноименных компонент калибровочных векторов. Предложена алгоритмическая процедура согласованной или так называемой многоскважинной (региональной) калибровки, ориентированная на выполнение этого условия. Рассмотрены возможности реализации этой процедуры в практических расчетах, предусматривающие, как автоматический, так и интерактивный режим ее проведения.
Послойная интерполяция калибровочных значений модельных параметров для определения компонент вектора (0)X должна учитывать особенности пространственного расположения Vk (T), калибровочных скважин k=0,1,2,…,n V. Как правило, (рис.7) эти скважины расположены в области (T) нерегулярно. В этой связи наиболее естественно осуществлять решение задачи интерполирования значений модельных параметров на триангуляционной сети, узлами которой являются калибровочные скважины.
В работе используется один из известных Рис. 7. Пример пространственного расположения калибровочных скважин. методов интерполяции на триангуляционных Скважины отмечены точками на карте сетях и предлагается его модификация, глубин залегания одной из формаций в позволяющая повысить устойчивость получаемых результатов на так называемых осадочном бассейне «узких» треугольниках.
(0) Значения компонент – вектора модельных параметров для прогнозной скважины V0 определяются в результате интерполяции с погрешностью. На основании оценки этой погрешности для каждого модельного параметра может быть задан интервал неопределенности. В совокупности эти интервалы образуют в пространстве модельных параметров многомерный параллелепипед 0 (0) неопределенности XX с центром в точке. В пространстве модельных данных P этот параллелепипед о тображается оператором прямой задачи в подмножество PP возможных модельных распределений эксцесса ГФД для скважины V0.
Дополнительный анализ этого подмножества позволяет, кроме прогноза P*0 (z), формируемого в скважине V0 соответствии с (24), дать оценку наихудшего и наилучшего вариантов возможного пространственного распределения эксцесса ГФД.
Под наилучшим распределением понимается нижняя оценка значений эксцесса ГФД в каждом узле сетки 0z, а под наихудшим – верхняя оценка этих значений.