CC BY
38
Вестник МГТУ, том 13, №4/1, 2010 г.
стр.793-809
УДК 519.87 : 551.1/.4
Эффективный подход к численному моделированию характеристик природных динамических систем
А.-В.И. Середа
Политехнический факультет МГТУ, кафедра высшей математики и программного обеспечения ЭВМ
Аннотация. Работа посвящена краткому изложению основных элементов методологии эффективного численного моделирования характеристик сложных природных динамических систем и рассмотрению некоторых практических аспектов ее реализации. Предлагаемый подход в концептуальном плане имеет достаточно общий характер. Вместе с тем, в данной работе он рассматривается для удобства изложения применительно к сформулированной для некоторого класса природных динамических систем задаче определения (прогнозирования) значений целевой характеристики исследуемой системы, которые эта характеристика принимает в настоящее время в заданной пространственной области. Методология прошла успешную апробацию на реальных задачах исследования геофлюидодинамических систем в осадочных бассейнах Земли.
Abstract. In the paper the methodology of effective numerical modeling of characteristics of natural dynamic systems has been discussed with reference to the problem of prediction of values of the investigated system target characteristic. This characteristic has the predicted values in the specified spatial area now. The methodology was successfully used for the problem of prediction of characteristics in sedimentary basins of the Earth.
Ключевые слова: природная динамическая система, эффективное моделирование, прямая и обратная задачи, численное моделирование, модельные параметры, анализ чувствительности
Key words: natural dynamic system, effective mathematical model, forward and inverse problem, numerical modeling, model parameters, sensitivity analysis
1. Введение
Под природной динамической системой Y(t) в работе понимается сплошная, пространственнонеоднородная по своим свойствам и составу, структурно определенная, в общем случае многокомпонентная материальная среда. Система Y(t) формируется во времени и заполняет собой некоторую ограниченную пространственную область f2(t), te [0,7], где t=0 - момент зарождения системы, t=T - настоящее время.
Развитие системы Y(t) во времени предполагает в общем случае возможность добавления в нее извне новых или исключения из системы уже входящих в нее структурных компонент, изменение свойств этих компонент вследствие протекающих в системе внутренних процессов различной природы (механических, физических, химических, термических и т.д.). В общем случае структура системы Y(t) и пространственные границы области £2(i) изменяются во времени.
Будем считать, что пространственная область £2(t) в любой момент времени te [0,7] может быть представлена как объединение пространственно упорядоченной совокупности конечного числа ограниченных, пространственно протяженных, односвязных и не имеющих общих точек подобластей (Oi(f)^Q(f), i=1,2,...,nn(t), nn(t)>1:
na(t)
£2(t) = urni(t) ; co1(t)naj(t) =0 V i±j, te [0,7]. (1)
При этом каждая из указанных подобластей заполнена физически однородной1 материальной средой. В дальнейшем будем называть эти подобласти слоями. Пространственная упорядоченность слоев mi(t), i=1,2,...,nn(t) определяет структуру исследуемой системы, а их количество и расположение в пространстве в общем случае зависят от времени. Систематизированное описание пространственно упорядоченной совокупности слоев сор), i=1,2,...,nn(t) будем называть в дальнейшем структурной моделью системы Y(t) и обозначать как Yg(t). В рамках сделанных предположений структурную модель Yg(t) можно называть слоистой. Для простоты будем предполагать, что слои юр), i=1,2,...,nn(t) не претерпевают структурных нарушений в области £2(t), а физические свойства материальной среды, как функции пространственных координат, в пределах слоя изменяются незначительно.
1 взаимная диффузия различных по характеристикам материальных сред предполагается пренебрежимо малой.
793
Середа А.-В.И. Эффективный подход к численному моделированию...
Допустим, что хронологически история формирования исследуемой системы на временном отрезке [0,7] от момента зарождения системы £(t) до сегодняшнего дня может быть задана с помощью временной сетки:
П = {T / 7 = T1.l+htl , htl>0, /'=1,2,...,wt; 7o = 0, Тщ = T}. (2)
Узлы сетки nt совпадают с моментами изменения внешних условий формирования системы, приводящими к структурным изменениям в ней. Эти структурные изменения заключаются либо в начале процесса вхождения в нее нового слоя, либо в исключении из нее некоторого числа уже имевшихся в системе слоев. На любом временном интервале 5Ti=(Ti.1,Ti), /=1,2,...,wt в системе может происходить (начиная с момента времени Ti-1 и заканчивая моментом времени Ti) процесс вхождения в систему только одного нового слоя. Для других принадлежащих системе слоев этот процесс считается завершенным. Исключение из системы некоторого конечного количества уже находящихся в системе слоев будем считать происходящим мгновенно. Свойства и пространственные границы всех слоев, образующих систему, могут изменяться во времени в силу протекающих в ней различных процессов.
Зададим для определенности в области Q(t) декартову2 систему координат OXYZ. Тогда каждой пространственной точке A(x,y,z)e £2(t), te [0,T] может быть поставлен в соответствие вектор a(x,y,z,t)e Es -вектор физических характеристик, определяющих локальные свойства системы X(t) в этой точке в момент времени t. Здесь Es - 5-мерное евклидово пространство, где s - количество указанных характеристик, часть из которых определяет свойства собственно материальной среды, образующей в совокупности исследуемую систему, а часть является следствием протекающих в системе процессов различного характера. В общем случае некоторая часть этих характеристик прямо или косвенно зависит друг от друга.
Состояние системы Х(0 в каждый момент времени te [0,T] определяется ее структурой и пространственными физическими полями - пространственными распределениями в области Q(t) количественных значений ее физических характеристик. В общем случае текущее состояние системы является результатом ее эволюционного развития в предшествующий промежуток времени. Предполагается, что в силу объективных причин, структура системы X(t) и значения ее физических характеристик потенциально доступны для прямых или косвенных измерений лишь в настоящий момент времени t=T.
Одной из задач при исследовании природных динамических систем является задача определения значений, которые принимает в настоящий момент времени интересующая исследователя физическая характеристика (целевая характеристика) системы X(t), в некоторой заданной подобласти области Д7). Такую задачу будем называть задачей прогнозирования характеристик природной динамической системы X(t). Термин "прогнозирование" используется здесь не во временном, а в пространственном смысле. Часто подобные задачи решаются посредством подходящей интерполяции или экстраполяции имеющихся данных о значениях целевой характеристики в требуемую пространственную область. Однако в ряде случаев такой подход к решению задачи неприемлем. В частности, это может иметь место в силу зависимости пространственных распределений физических характеристик в области EMT) не только от конкретного состояния системы, но и от истории ее предшествовавшего развития. В этом случае задача должна решаться в контексте эволюционного развития системы от момента ее зарождения до настоящего времени, что является существенным отличием от традиционных задач обработки и интерпретации данных наблюдений.
2. Постановка задачи прогнозирования характеристик природной динамической системы
Будем считать справедливыми следующие утверждения:
а) Скорость формирования системы X(t) и скорости протекающих в ней при этом процессов очень малы по сравнению с процессами, проходящими в масштабе реального времени. Как следствие, временной промежуток [0,T] слишком велик для того, чтобы имелась возможность непосредственно наблюдать весь жизненный цикл системы от ее зарождения до настоящего времени.
б) Система X(t) потенциально доступна для непосредственного комплексного изучения (посредством прямых или косвенных измерений) лишь в настоящий момент времени3 (при t=T).
в) Значения компонент вектора a(x,y,z,T) - физических характеристик исследуемой системы, как правило, недоступны для прямых измерений во всей области £2(T). Такие измерения возможны лишь для отдельных характеристик и лишь для, вообще говоря, конечного множества точек A(x,y,z)e £2(T).
г) Существуют комплексные косвенные методы исследования, позволяющие с определенной степенью точности определить временной отрезок [0,T] "жизни" системы, задать временную сетку nt, позволяющую хронологически отобразить историю формирования системы X(t). В результате
2 Тип используемой системы координат зависит от специфики пространственного строения исследуемой системы или свойств изучаемых процессов.
3 Впрочем, именно этот момент времени и имеет обычно наибольшее практическое значение.
794
Вестник МГТУ, том 13, №4/1, 2010 г.
стр.793-809
оказывается возможным построить структурную модель системы Xg(0 в любой момент времени te [0,7] и определить для каждого из временных интервалов ST1=(T1.1,T1), i=1,2,...,nt возможные (стартовые) свойства слоя (d(t), i=l,2,...,na(t) при вхождении его в систему.
д) Процессы, протекающие в системе, в достаточной мере изучены, что позволяет при необходимости использовать ряд упрощенных эмпирических законов, количественно описывающих изменение ее свойств с течением времени.
Для определенности будем считать, не умаляя общности, что пространственная область £2(t) здесь и далее задается в каждый момент времени te [0,7] (в выбранной декартовой системе координат4) условиями:
xmin—xmax ,
| ymm—x— ymax, (3)
l 0 — Z— Zmax(x, y, t).
В качестве основной задачи исследования рассмотрим следующую проблему.
Предположим, что структурная модель исследуемой системы Xg(t) известна для любого te [0,7].
Пусть имеется конечное множество пространственно локализованных подобластей Vk<^Q(T), £=0,1,2,...,nV, каждая из которых пересекается с одним и тем же множеством слоев структурной модели Yjg(T). Для определенности будем считать, что подобласти Vk представляют собой вертикальные одномерные разрезы области Q(T). Тогда
Vk = {(Xk,yk,,z)eQ(T) / 0 — z— ZmaxfeybT)}, £=0,1,2,...,nv.
Для каждого вертикального разреза Vk<^Q(T), £=1,2,...,nv определена одномерная сетка:
%z = {zki / Zki = Zk(i-1)+hzki , hzki>0, i=1,2,.,nzk; Zk0 = ° Zknzk — Zmax(VЬT)},
в узлах Zki которой известны (в результате прямых или косвенных измерений) значения P (xk,yk,zki,T) -целевой характеристики в настоящее время. Эти значения будем в дальнейшем называть наблюдаемыми или полевыми данными, а совокупность значений P (xk,yk,zki,T), i=1,2,.,nzk обозначать как Pk(z) и называть дискретным распределением целевой характеристики для вертикального разреза Vk.
Вертикальный разрез V0<^Q(T) будем называть прогнозным. Для этого вертикального разреза полевые данные не заданы.
Сформулируем теперь основную задачу.
Основная задача
Пусть для каждого вертикального разреза Vk<^Q(T) известно полевое дискретное распределение целевой характеристики P k(z), k=1,2,.,nv.
Требуется для заданного (прогнозного) вертикального разреза V0<^Q(T) определить дискретное распределение P0(Z) в узлах сетки
n0,z = {z0i / z0i = Z0i-1+hz0i , hz0i>0, ^1,2,‘ • *,nz0; z00 = ° Z0nz0 — zmax(V0,T)}.
При решении Основной задачи будем учитывать следующие предположения.
1. Получение требуемого результата посредством интерполяции известных полевых распределений Pk(z), k=1,2,...,nv целевой характеристики в прогнозный вертикальный разрез не представляется возможным, поскольку пространственное распределение этой характеристики в области Q(T) является результатом эволюционного развития системы X(t) с момента ее зарождения до настоящего времени. Решение задачи должно происходить в контексте всего периода эволюционного развития системы.
2. В общем случае практически отсутствует достаточно точная априорная количественная информация о свойствах системы в настоящее время, количество nv вертикальных разрезов Vk с полевыми данными о дискретных распределениях Pk(z) обычно относительно невелико, а представленные в них измерения значений характеристики, как правило, неравноточные и часто имеют большую погрешность.
3. Может быть построена (например, на основе фундаментальных законов сохранения) физически интерпретируемая математическая модель, позволяющая находить (с учетом замечания 1) P(x,y,z,T) -модельное распределение целевой характеристики в области £2(T). Как правило, практическое построение модельного пространственного распределения P(x,y,z,T) осуществляется в результате аппроксимации такой модели дискретной моделью (например, разностная аппроксимация краевой задачи
4 Ось OZ считается направленной вниз.
795
Середа А.-В.И. Эффективный подход к численному моделированию...
для уравнения математической физики) с ее последующим численным исследованием на ЭВМ. Представим итоговую зависимость P(x,y,z,T) от модельных параметров в явном виде:
Px = GX), X^X, PxeP, (4)
где Px= PX(x,y,z,T)eP<^En - модельное распределение целевой характеристики; £eXcEm - вектор варьируемых модельных параметров; G(x) - оператор, осуществляющий отображение из X -пространства модельных параметров в P - пространство модельных распределений.
Здесь и далее к варьируемым модельным параметрам необходимо отнести некоторую часть таких свойств среды, которые, с одной стороны, существенны для моделируемого процесса, а с другой стороны, количественная оценка которых затруднительна (известен, например, лишь диапазон их возможных значений).
Задачу (4) будем в дальнейшем называть задачей прямого моделирования или прямой задачей, а оператор G(X) - оператором прямой задачи.
4. Непосредственное получение P0(z) - требуемого прогноза дискретного распределения целевой характеристики для заданного вертикального разреза V0(x0,y0,z)^Q(T) с помощью операторного соотношения (4) практически неосуществимо по нескольким причинам.
Во-первых, даже если предположить, что построение точной модели изучаемых процессов возможно, эта модель, как правило, оказывается для реальных природных систем чрезвычайно громоздкой и сложной не только для аналитического, но и для численного исследования.
Во-вторых, из-за очевидного недостатка сведений о фактических свойствах среды объективно исключается возможность получения достаточно точных результатов моделирования изучаемых процессов даже при наличии точной модели и приемлемых методов ее исследования.
5. В пределах каждого слоя ф(t), /=1,2,...,wn(t) структурной модели Xg(0 исследуемой системы для значений модельных параметров, как правило, характерна незначительная вариабельность. При этом в целом изменение этих значений в пространстве и во времени закономерно с физической точки зрения.
3. Принципы, положенные в основу методологии численного моделирования
При разработке методологии численного решения Основной задачи предлагается исходить из следующих базовых принципов.
Принцип эволюционного развития
Постановка и решение прямой задачи (4) должны осуществляться в контексте эволюционного развития системы Х(0 с учетом динамики изменения ее структуры и свойств образующей ее материальной среды.
Принцип эффективного моделирования
Для моделирования исследуемого процесса должны использоваться эффективные модели. Под эффективной моделью понимается такая физически обоснованная модель исследуемого процесса, построение которой осуществляется для обеспечения надежного моделирования заданной (целевой) характеристики процесса. При этом перед эффективной моделью не ставится задача точного и адекватного описания моделируемого процесса в целом. Надежность моделирования понимается в том смысле, что эффективная модель должна быть потенциально настраиваемой (калибруемой) на имеющиеся полевые данные о значениях целевой характеристики. Другими словами, каковы бы ни были полевые данные, всегда возможно выбрать такие значения модельных параметров (из множества их возможных значений), при которых моделируемые значения целевой характеристики (модельные данные) будут в достаточной степени близки к полевым данным. По существу эффективная модель становится функционально адекватной лишь после ее калибровки (параметрической идентификации модели).
Введенный принцип эффективного моделирования носит общий характер и справедлив не только по отношению к математическим моделям, но и к любым другим видам моделей. В этой связи необходимо особо выделить эффективные численные модели. Под численной моделью, построение которой является необходимым при реализации численного моделирования, понимается совокупность вычислительных операций, осуществляемых над исходными данными, являющаяся алгоритмическим результатом применения выбранного численного метода анализа исследуемой модели. Эффективная численная модель понимается в смысле общего определения эффективной модели. Численное моделирование с использованием эффективных численных моделей будем называть эффективным. Большое значение в качественной оценке эффективного численного моделирования имеет устойчивость получаемых результатов. Важно отметить, что собственно погрешность численного анализа, которая складывается из погрешности модели, метода и вычислительной
796
Вестник МГТУ, том 13, №4/1, 2010 г.
стр.793-809
погрешности, не играет при этом существенной роли. Необходимая точность (при обеспечении устойчивости получаемых в результате эффективного численного моделирования результатов) всегда может быть надежно обеспечена в результате калибровки эффективной численной модели.
По сравнению с точными моделями, эффективные модели, как правило, существенно более просты с точки зрения возможности их теоретического и численного анализа, что является их основным и несомненным преимуществом. С другой стороны, они имеют целенаправленный характер и, в общем случае, не могут служить основанием для физически адекватного изучения моделируемых процессов в целом.
В качестве существенного условия для эффективных моделей целесообразно принять обеспечение возможности естественной физической интерпретации в терминах предметной области как самих модельных описаний, так и их (модельных) параметров.
Необходимость осуществления параметрической идентификации (калибровки) эффективной модели предполагает постановку и решение обратной задачи, которая в общем случае может быть сформулирована, например, в следующем виде:
Задача х
Пусть в области £2(T) известно полевое распределение P (x,y,z,T)eP целевой характеристики. Требуется найти вектор х eX, для которого выполнено условие:
G(X) = P ,X eX, P*eP. (5)
Обратная задача (5), как правило, не является корректной, и для ее решения требуется применение тех или иных регуляризирующих подходов.
Кроме того, необходимо обеспечить приемлемую вычислительную трудоемкость решения задачи (5). В этой связи естественным является следующий принцип, принятый в данной работе.
Принцип минимальной сложности
В соответствии с этим принципом используемые эффективные модели должны быть предельно просты с вычислительной точки зрения. В связи с этим из принятого класса эффективных моделей необходимо выбрать модель по возможности минимальной вычислительной сложности. Реализацию такой модельной идентификации возможно осуществить посредством рациональной минимизации количества варьируемых модельных параметров прямой задачи за счет исключения из их числа тех параметров, влияние изменения значений которых на результат моделирования незначительно. С этой же целью необходимо осуществить построение структурной модели Xg(0, te [0,T] минимально допустимой сложности, что достигается за счет рационального уменьшения nn(t) - общего числа входящих в нее структурных элементов - слоев (d(t), /-1,2,.. ,,nn(t).
Наконец, следует отметить еще один существенный принцип, используемый при построении предлагаемой в работе методологии численного прогнозирования.
Принцип пространственной локализации
Физические свойства системы, ее структурная организация могут объективно способствовать специфической пространственной ориентации или локализации протекания в ней исследуемого процесса. В этом случае и сбор данных, и моделирование желательно осуществлять в соответствии с имеющейся пространственной спецификой изучаемого процесса.
Кроме того, по объективным причинам полевые данные (безотносительно к особенностям исследуемых процессов) могут быть пространственно локализованы или пространственно ориентированы. И в этой ситуации моделирование должно, по возможности, осуществляться с учетом пространственной специфики имеющихся полевых данных.
Так, в данной работе (при формулировке Основной задачи) было сделано предположение, что полевые данные о значениях моделируемой характеристики представлены только в пространственно локализованных вертикальных разрезах - подобластях Fkofl(T), k=0,1,2,...,nV. То есть соответствующие полевые данные о дискретных распределениях Pk(z), £=1,2,.,nV пространственно ориентированы вдоль оси OZ в выбранной декартовой системе координат.
В этой связи часто оказывается целесообразным произвести декомпозицию исходной математической модели на так называемую одномерную математическую модель, ориентированную на моделирование "проекции" исследуемого процесса на ось OZ и двумерную - ориентированную на моделирование "проекции" процесса на горизонтальную плоскость OXY. Одномерную модель, полученную в результате такой декомпозиции, будем называть вертикальной моделью, а двумерную - латеральной. При этом неучет в вертикальной модели латеральных составляющих процесса компенсируется за счет введения в модель дополнительных эффективных модельных параметров. Аналогичным образом осуществляется эффективный учет вертикальных составляющих процесса в латеральной модели.
797
Середа А.-В.И. Эффективный подход к численному моделированию...
Следует отметить, что в ряде случаев оказывается возможным ограничиться анализом лишь вертикальной модели применительно к каждому вертикальному разрезу, что и будет предполагаться в дальнейшем.
Для вертикальной модели прямую задачу (4) можно переформулировать в виде:
Pck (z)= GZ(X\ X^eXz, PzkePz, k=1,2,...,«v, (6)
где Gz(x) - оператор, осуществляющий отображение из Xz - пространства модельных параметров вертикальной задачи в Pz - пространство модельных распределений для вертикальной модели. Pxk(z)=Pxk(xk,yk,z,T)ePz - модельное распределение целевой характеристики, а - вектор
варьируемых модельных параметров для вертикального разреза Vk<^Q(T), k=1,2,...,nv.
Задачу (6) будем в дальнейшем называть одномерной задачей прямого моделирования или вертикальной прямой задачей, а оператор Gz(%) - оператором вертикальной прямой задачи.
Соответственно, калибровку вертикальной модели теперь естественно проводить отдельно для каждого вертикального разреза Vk^Q(T), k=1,2,...,nv, решая обратную задачу в следующей постановке:
Задача у,
Пусть для вертикального разреза Vk<^Q(T) известно полевое дискретное распределение Pk(z)ePz, целевой характеристики. Требуется найти такой вектор X'k)*£ Xz, для которого выполнено условие:
Gz(Xk)*) = P*k(z), fk)*e-Xz, P*k ePz. (7)
В отношении задачи (7), с точки зрения ее некорректности, могут быть сделаны те же замечания, что и для задачи (5). Однако в вычислительном отношении нахождение численного решения задачи (7), безусловно, будет существенно менее трудоемким. Впрочем, как и задачи (6) по сравнению с задачей (4).
В контексте задачи (7) вертикальные разрезы Vk(xk,yk,z)c^(T), k=1,2,...,wv будем называть в дальнейшем калибровочными разрезами, а вектора х®* - соответствующими калибровочными векторами модельных параметров.
Отметим в заключение, что рассмотренные выше принципы эволюционного развития, эффективного моделирования, минимальной сложности и пространственной локализации характеризуют в совокупности концепцию эффективного моделирования, которая положена в основу предлагаемого эффективного подхода к численному моделированию характеристик природной динамической системы.
4. Калибровка модели. Регуляризирующий подход к решению обратной задачи
Как уже отмечалось, сформулированные выше задачи о калибровке (5) и (7) являются в математическом отношении обратными по отношению к задачам (4) и (6) соответственно.
Подобного рода обратные задачи часто возникают в естественнонаучных исследованиях. Они обладают рядом принципиальных математических особенностей, которые необходимо учитывать при разработке или выборе методики их численного анализа. Некоторые важные для дальнейшего положения и свойства обратных задач рассматриваются ниже.
Определение корректной задачи, как известно, ввел в начале 30-х годов прошлого века Ж. Адамар (Hadamard, 1932). Запишем это определение в принятых выше обозначениях.
Определение 1
Задача определениях, удовлетворяющего условию:
G(x) = P‘, X вХ, P'e P, (8)
где X и P - некоторые заданные метрические пространства, является корректно поставленной (корректной) по Адамару, если для нее выполнены следующие три условия:
1. задача (8) имеет решение на всем пространстве P;
2. решение задачи (8) единственно;
3. решение задачи (8) устойчиво.
Последнее условие означает, что малым возмущениям в правой части уравнения (8) будут соответствовать малые возмущения в решении задачи. Очевидно, что строго говоря задачи (5) и (7) в общем случае некорректны. Причем без дополнительных предположений для них не выполнено ни одно из условий определения (8). Действительно, рассматривая задачу (5) (для второй задачи рассуждения аналогичны), можно отметить следующее:
- В общем случае нет возможности гарантировать существование решения задачи на всем пространстве P. Безотносительно к конкретным математическим свойствам оператора G(x) одной из причин этого для конкретных природных динамических систем является, например, обязательная
798
Вестник МГТУ, том 13, №4/1, 2010 г.
стр.793-809
физическая интерпретируемость результатов моделирования с помощью оператора прямой задачи G(%), что принципиально исключает существование решений задачи (5) для тех Р е Р, физически оправданная интерпретация которых невозможна.
- Если решение задачи (5) существует, то его неединственность в общем случае объективно обусловлена природой исследуемых процессов, когда различные физические условия (определяемые различными значениями модельных параметров) могут приводить к одному и тому же или близкому конечному результату.
- Наконец, не может быть гарантирована и устойчивость возможных решений задачи (5), что объясняется наличием в общем случае среди модельных параметров таких, любое допустимое изменение значений которых оказывает незначительное воздействие на изменение решений прямой задачи. Будем называть такие параметры малочувствительными. Очевидно, что чувствительность к изменению значений модельных параметров является локальным свойством оператора прямой задачи G(x). Поэтому необходимость отнесения того или иного модельного параметра к множеству малочувствительных параметров и общее количество таких параметров в различных точках пространства модельных параметров могут изменяться.
Таким образом, обратные задачи (5) и (7) в общем случае следует отнести к классу некорректных задач. В ходе дальнейших рассмотрений будем следовать в основном работам (Гончарский, 1987; Бакушинский, Гончарский, 1989).
Сформулируем понятие условно корректной задачи, введенное А.Н. Тихоновым (1943): Определение 2
Задача (8) является условно корректной при выполнении следующих условий:
1. Известно, что решение задачи (8) существует и принадлежит ХсХ - некоторому известному подмножеству Х;
2. Оператор G осуществляет взаимно однозначное отображение Х0 на Р0сР;
3. Оператор G"1 непрерывен на Р0.
Таким образом, для условно корректной задачи не требуется разрешимость во всем пространстве Р, что имеет место для задачи, корректной по Адамару. Требование непрерывности G"1 на пространстве Р заменяется требованием непрерывности G"1 на его подмножестве Р0 - образе подмножества Х°с Х.
Множество Х в определении условно корректной задачи называют множеством корректности. В частности, в качестве множества корректности в задаче (8) может быть выбрано (Гончарский, 1987) любое компактное подмножество Х.
Формирование множества корректности Х и рассмотрение задачи (8) в классе условно корректных задач с множеством корректности Х0
является ключевой проблемой в вопросе построения метода решения этой задачи. В частности, для решения некорректных задач на компактных множествах корректности в работах В.К. Иванова (1963; 1970) разработан метод так называемых квазирешений, согласно которому квазирешением задачи (8) на множестве корректности Х называется элемент хе еХ, минимизирующий на этом множестве значение IIРх - Ре* II . То есть
IIPxs* - Ре* II = inf II Рх - Ре* II. (9)
*еХ 0
Здесь под Ре* будем понимать заданное распределение целевой характеристики, а индекс е означает, что это распределение задается с некоторой погрешностью.
Нетрудно показать, что алгоритм построения квазирешения X - приближенного решения задачи (8) в соответствии с (9) - является регуляризирующим алгоритмом. Его существенной особенностью является то, что построение осуществляется только с использованием Ре* и не требует знания е.
Возможен и несколько модифицированный (Гончарский, 1987) подход к построению
квазирешения условно корректной задачи на компактном множестве корректности, который в большей степени соответствует практическим требованиям, предъявляемым к решениям обратной задачи в данной работе. Этот подход заключается в том, что в отличие от (9), в качестве хе еХ0 - приближенного решения задачи (8) выбирается любой вектор, для которого выполнено условие:
I I РХе. - Ре* I I < е . (10)
Множество {%еХ/ I I РХе* - Ре I I < е} не пусто, поскольку оно всегда включает в себя, по крайней мере, одну точку - точное решение задачи (8). Поэтому удовлетворяющее (10) приближенное решение Хе еХ0 всегда может быть построено. При этом, очевидно, что метод (10) построения хе также является регуляризирующим алгоритмом, так как выполнено условие:
799
Середа А.-В.И. Эффективный подход к численному моделированию...
lim Xe* = X-
£——0
В методе (10) требуется априорное знание не только Pe, но и £. Вместе с тем для метода (10) практически снимается проблема существования и единственности квазирешения. Решения задачи (8), получаемые в соответствии с методом (10), будем (Бакушинский, Гончарский, 1989) называть в дальнейшем £-квазирешениями.
Существенно отметить также, что для случая конечномерных пространств X и P и ограниченного множества ХсХ методы (9) и (10) могут быть реализованы с помощью известных методов минимизации функционалов. В том числе, в ряде случаев (Бакушинский, Гончарский, 1989) - в рамках стандартной схемы метода наименьших квадратов.
В результате можно сделать следующие практически важные выводы:
- задачи (5) и (7) в общем случае являются некорректными задачами;
- задачи (5) и (7) могут быть отнесены к условно корректным задачам с соответствующим множеством корректности X0,
если обеспечить выполнение условий 2 и 3 Определения 2;
- для решения задач (5) и (7), сформулированных в классе условно корректных задач, в качестве регуляризирующего метода могут быть использованы методы (9) или (10).
С учетом проведенных рассмотрений скорректируем формулировки прямой задачи (4) и задачи о калибровке (5). При этом будет также учитываться наличие погрешностей в полевых данных и фактически имеющие место ограничения на выбор вектора X еХ1.
Введем в рассмотрение множество XсХ, которое задает множество векторов X со значениями модельных параметров из некоторой естественной области их возможных значений, определяемой, например, из физических предположений, и перепишем прямую задачу (4) в форме:
Px = G(X), XeX°cX, PxgPVP, а!
где P.X=PX(x,y,z,T)eP°<zP<zEn - модельное распределение целевой характеристики; %^X><zX<zEm - вектор варьируемых модельных параметров; G(x) - в общем случае нелинейный оператор прямой задачи, однозначно отображающий X0 - замкнутое подмножество конечномерного евклидового пространства модельных параметров X, на P0 - замкнутое подмножество конечномерного евклидового пространства модельных распределений P. Компонентами вектора модельных параметров X являются значения существенных для моделируемого процесса характеристик физических свойств среды.
Множество ХоУ может быть задано, в частности, как многомерный "параллелепипед" в векторном пространстве X с помощью указания диапазонов возможного изменения значений каждого из модельных параметров:
Ximm < Xi < Ximax, i=1,2,...,m. (12)
Будем предполагать в дальнейшем непрерывность оператора G(x) на X. С физической точки зрения, непрерывность оператора G(x) означает, что небольшим вариациям в значениях компонент вектора модельных параметров Xе X соответствуют небольшие же изменения компонент модельного распределения PxeP°. В общем случае, обеспечение непрерывности оператора G(x) может привести к необходимости введения дополнительных ограничений на выбор множества X <zX .
Пусть теперь, как и раньше, P =P (x,y,z,T) - полевые данные о пространственном распределении целевой характеристики в области £2(T), полученные в результате прямых или косвенных измерений. Это распределение будем называть наблюдаемым, или полевым, и обозначать для краткости P.
Будем считать, что измерения проводятся в узлах некоторой пространственной трехмерной прямоугольной в общем случае неравномерной сетки ®xyz(T):
®xyz(T) = {(xbУj, zk / xi=xi-1+hxb hxi>0, i=1,2,.,^x; x0 — xmin; xnx < xmax;
yj =yj-1 +hyj, hyj>0, j=1,2,.,ny; y0 — ymin; yny < ymax; 0-3)
zk=zk -1 ^^z^ hzk >0, &=1,2,. • • ,nzxy; z0 = ^ znzxy < zmax(xk,yb T)}.
В этом случае не умаляя общности можно считать, что P представляет собой конечномерный вектор. Для обеспечения сопоставимости данных будем предполагать, что Px - модельное распределение целевой характеристики, получаемое из решения прямой задачи (11), также представлено в узлах сетки fflxyz(T) для VXeX0, что всегда может быть обеспечено за счет соответствующей организации вычислительного процесса.
Очевидно, что для произвольного XeX модельные и полевые данные распределения целевой характеристики в общем случае не будут совпадать с достаточной степенью точности в метрике пространства P. Вместе с тем прямая задача (11) лишь тогда является приемлемой для того или иного
800
Вестник МГТУ, том 13, №4/1, 2010 г.
стр.793-809
заданного вертикального (глубинного) разреза, когда указанное совпадение имеет место. Достичь этого в рамках заданного математического описания возможно лишь варьируя значения модельных параметров. Как следствие, возникает так называемая задача о калибровке модели. Здесь может быть использован также термин (Алифанов и др., 1988) "параметрическая идентификация" модели. Указанная задача заключается в определении такого вектора модельных параметров х*еХ, для которого полученное в соответствии с (11) модельное распределение целевой характеристики Px будет с достаточной степенью точности совпадать с известным полевым распределением P. Полевые данные для простоты будем считать равноточными. Введем в рассмотрение векторную функцию р(х, P):
И(Х,Р) = Px - P. (!)
Тогда задачу о калибровке модели (11) можно сформулировать как задачу нахождения вектора х'еХ, для которого выполнено условие:
llp(z‘,P‘) II < s, (15)
где S > 0 - требуемая точность согласования Px и P .
Условие (15) является принципиальным, поскольку как модельные, так и полевые данные могут быть определены лишь с некоторой погрешностью. Следовательно, абсолютно нецелесообразно требовать их абсолютной согласованности.
В дальнейшем будем считать, что прямая задача (11) вполне адекватна исследуемому процессу, понимая под этим, что вектор х еХ, для которого выполнено условие (15), существует.
Введем в рассмотрение US(P0) - S-расширение множества P0:
US(P0) = (PeP / min I I цХ P)| I < S}. (16)
XeX0
Очевидно, что для обеспечения адекватности прямой задачи в указанном выше смысле достаточно потребовать, чтобы вектор полевых данных P принадлежал S-расширению множества P0:
Pe U3(P°). (17)
Для удобства дальнейшего рассмотрения переформулируем теперь задачу (5) с учетом принятых выше допущений и обозначений и формулировки задачи (11).
Задача о калибровке
Пусть имеется численная модель целевой характеристики, задаваемая операторным уравнением
Px=G(x), XeX0, PxeP0, (18)
где ХеХ° - вектор модельных параметров; PxeP0 - вектор расчетных (модельных) значений целевой характеристики; G(x) - в общем случае нелинейный непрерывный на X оператор, взаимно однозначно отображающий X в P0; X0<zX - замкнутое, ограниченное подмножество конечномерного евклидового векторного пространства X; P0cP - замкнутое, ограниченное подмножество конечномерного евклидового векторного пространства P.
Пусть далее PE*e Ug(P°) - вектор наблюдаемых значений целевой характеристики (полевые данные). Будем считать измерения равноточными и оценивать погрешность задания PE значением е>0, понимая под этим выполнение условия:
I I P/- P* I I < е, (19)
где P - вектор с точными значениями полевых данных.
Требуется определить Хв eX - вектор модельных параметров, удовлетворяющий условию:
I I МХв>Х I I < S , (20)
где векторная функция ^(х) определяется в соответствии с соотношением
M/,P/) = Px - Pe*. (21)
Значение S>0 определяется с учетом е - погрешности полевых данных и, возможно, вычислительной погрешности или задается экспертно.
Очевидно, что формулировка задачи (18-21) предполагает построение S-квазирешений задачи, обратной к задаче (11). В качестве метода ее решения может быть в принципе использован любой подходящий метод минимизации функции рассогласования:
F(X) = I I MX,Pe*) I I 2 (22)
801
Середа А.-В.И. Эффективный подход к численному моделированию...
на множестве X0 При этом работа метода может быть завершена, как только будет выполнено условие (20).
Ключевой проблемой в вопросе построения метода решения задачи (18-21) является формулировка ее в классе условно корректных задач с множеством корректности X0. Построение такого множества может быть осуществлено в соответствии с условиями (12) при исследовании конкретной системы исходя из общих свойств и специфических особенностей конкретного оператора G(%). Как правило, вначале множество X0
определяется априорно на основании имеющейся информации или предположений о возможных диапазонах изменения значений модельных параметров. Затем, по мере необходимости, оно уточняется алгоритмически в процессе решения задачи (18-21). Так, например, при построении множества X0
может оказаться полезным сужать диапазоны возможного изменения значений (вплоть до их фиксации) тех модельных параметров, к изменению значений которых модель малочувствительна, что будет способствовать повышению устойчивости получаемых решений обратной задачи.
Вместе с тем, следует отметить, что выполнение условия (20), конечно, не обеспечивает в общем случае единственности полученного решения формулируемой задачи о калибровке. Это условие позволяет лишь определить одно из ее 5-квазирешений, принадлежащих некоторому подмножеству множества X0,
для всех точек которого выполнено условие (20). Кроме того, в случае неточного построения множества корректности X0 (например, G(x) осуществляет однозначное отображение X0 в P0, а обратный оператор неоднозначен) таких подмножеств множества X0 может быть больше, чем одно. В этом случае выбор наиболее подходящего из них в смысле того или иного критерия должен обеспечиваться, например, заданием дополнительной априорной информации, уточняющей множество X0, обеспечивая исключение из него всех указанных подмножеств, кроме требуемого.
5. Снижение размерности пространства модельных параметров
В соответствии с принципом эффективного моделирования, упрощение модели не должно приводить к нарушению требований физической интерпретируемости используемых построений и возможности обеспечения требуемого уровня согласованности модельных и полевых данных. В этой связи существенным фактором является определенная сбалансированность сложности исследуемых природных процессов, сложности и детальности модельных описаний, требуемой точности модельных данных, точности, количества и структуры полевых данных. Отсутствие совокупной сбалансированности указанных элементов может приводить к ухудшению общих свойств получаемых модельных построений.
Так, формальная избыточность полевых данных может способствовать неоправданному усложнению модельных описаний. При этом качественного улучшения результатов моделирования не происходит при явном увеличении трудоемкости вычислительного процесса. В этом случае предварительный анализ и дифференциация полевых данных могут существенно снизить их формальную избыточность и повысить качество моделирования.
С другой стороны, сложные модельные построения при очевидной недостаточности или большой погрешности полевых данных будут порождать множественность трудно интерпретируемых и далеких от действительности результатов моделирования.
Общая оценка сложности модели предполагает чаще всего совокупный анализ сложности ее теоретического и численного исследования. Главной целью такого анализа в данной работе является обеспечение фактической эффективности численного исследования. Поэтому в качестве оценки сложности используемых моделей примем трудоемкость их численного анализа. Это оправданно с точки зрения необходимости разработки методологии численного прогнозирования, обеспечивающей возможность решения поставленной задачи в приемлемое время. Для сложных природных систем достижение этой цели является непростой и актуальной проблемой.
Таким образом, основной целью построения моделей минимальной сложности в данной работе будем считать максимально возможное снижение трудоемкости их численного анализа при обеспечении физической интерпретируемости используемых моделей и требуемого уровня согласованности модельных и полевых данных. Последнее будем использовать в качестве критерия сбалансированности сложности модельных описаний и полевых данных. При этом оказывается возможным построение процедур (Мадатов, Середа, 2003) поэтапного упрощения модели до тех пор, пока обеспечено выполнение указанного критерия.
В этом контексте упрощение численных моделей может осуществляться главным образом в трех направлениях:
- за счет выбора модельных описаний предельно допустимой простоты;
- за счет выбора экономичных численных методов;
- за счет снижения размерности пространства модельных параметров в результате уменьшения их количества для выбранного класса моделей.
802
Вестник МГТУ, том 13, №4/1, 2010 г.
стр.793-809
Проблема выбора модельных описаний предельной простоты не может быть конструктивно осуществлена в общем случае. Для этого как минимум необходимо иметь информацию о конкретных свойствах и структуре исследуемой системы, о типах и свойствах протекающих в ней процессов. Аналогично, выбор экономичных численных методов может быть осуществлен лишь на основании известных свойств и особенностей моделей, используемых для исследования той или иной конкретной природной системы.
Рассмотрим здесь более подробно третье направление - уменьшения количества модельных параметров в выбранном классе моделей. Принципиальное обсуждение этого направления может быть проведено безотносительно к особенностям конкретной природной системы и используемых для ее исследования моделей.
Используя принятые ранее обозначения, предположим, что для исследуемой в пространственной области Q(t) природной динамической системы £(t) построена структурная модель £g(t), представляющая собой упорядоченную в пространстве совокупность попарно соприкасающихся слоев (01(t)cQ(t), i=1,2,...,nn(t), na(t)>\, te[0,T]. Пусть, кроме того, задано конечное множество варьируемых модельных параметров и построен оператор прямого моделирования G(%), хеХ<zX<zEm. Каждому слою (01(T)cQ(T), i=1,2,...,na(T), na(T)>1, входящему в структурную модель £g(T), соответствует конечное число модельных параметров - компонент вектора %eX0. Для простоты примем, что каждому слою соответствует тх модельных параметров. Тогда т - общее число модельных параметров, будет определяться формулой:
т = mx-na(T).
Будем считать далее, что Х°сХ - множество корректности для оператора G(%). Пусть также задан параметр с>>0, определяющий требуемую точность согласования модельных и полевых данных, и известен вектор полевых данных РЕ*е Ug(P0)cPcEn. Функционал рассогласования F(%) определим в соответствии с формулой (22).
Предположим, что в указанных условиях может быть обеспечено необходимое согласование модельных Рх и полевых данных РЕ . То есть в результате решения задачи о калибровке модели в постановке (18-21) может быть определен вектор £Е еХ0, такой, что выполнено условие:
F(x:)<S2. (23)
Это означает, что потенциально используемая модель может быть упрощена, поскольку имеющие место модельные построения возможно избыточны с точки зрения принятого критерия согласованности реультатов моделирования с полевыми данными.
При построении процедуры упрощения численной модели будем исходить из того, что модель оптимальна по сложности, если, при прочих равных условиях, она обеспечивает достижение заданной точности калибровки (значение функционала рассогласования) при минимальной размерности пространства модельных параметров. Иными словами, задача заключается в определении минимального количества варьируемых модельных параметров, достаточных для обеспечения требуемого уровня согласованности модельных и полевых данных.
Рассмотрим схему возможной двухэтапной стратегии уменьшения размерности пространства модельных параметров.
На первом этапе сокращение размерности пространства модельных параметров осуществляется в результате эффективного упрощения структурной модели £g(T) исследуемой системы. Основная идея такого упрощения заключается в пространственном объединении5 двух входящих в структурную модель системы и пространственно соседствующих слоев (01(T) и (Ol+1(T) (1< i < nn(T)-l) в один слой. В результате такого объединения общее количество элементов структурной модели Xg(T) уменьшается на единицу. Соответственно, количество модельных параметров снижается на тх единиц. Для скорректированной соответствующим образом прямой задачи решается задача о калибровке (18-21). Если в результате решения этой задачи будет определен вектор £Е еХ0, такой, что выполнено условие (23), то объединение слоев (0,(T) и щ1+1(Г) считается обоснованным и не отменяется. Если требуемый уровень согласованности модельных и полевых данных не был достигнут, то объединение слоев (01(T) и (01+1(T) считается необоснованным и должно быть отменено. Рассмотренная процедура проверки обоснованности попарного объединения слоев последовательно осуществляется для всех 1< i < nn(T)-1. Будем называть это циклом процедуры. Циклы повторяются до тех пор, пока количество слоев в структурной модели £g(T) системы после очередного обоснованного объединения не окажется равным единице, либо после завершения выполнения цикла, в котором ни одно объединение не оказалось 5
5 Если такое объединение допустимо с физической точки зрения.
803
Середа А.-В.И. Эффективный подход к численному моделированию...
обоснованным. Если в результате количество слоев структурной модели системы уменьшится на n1, то общее количество модельных параметров снизится на величину mx-n1.
Второй этап предлагаемой стратегии уменьшения реальной размерности пространства модельных параметров заключается в определении и исключении из множества варьируемых так называемых малочувствительных параметров.
Количественная оценка Sj чувствительности6 модельного параметра Xj , j=\,2,...,m в точке ХеХ может быть определена, например, на основании относительного изменения функции рассогласования F(X) при некотором (нормированном) приращении значения
S = (F(X+&) - F(X))/F(X),
где e1 - единичный вектор, все компоненты которого равны нулю кроме j-й компоненты, которая равна единице; 8>0 - приращение значения параметра Х\-
Выделенным малочувствительным параметрам присваиваются наиболее вероятные значения, и в последующем при калибровке модели они учитываются как постоянные величины. Как уже отмечалось, чувствительность модели к изменению значений модельных параметров является ее локальным свойством в пространстве модельных параметров. Поэтому разбиение множества модельных параметров на подмножества чувствительных и малочувствительных параметров остается в общем случае справедливым лишь в некоторых малых окрестностях точек, в которых это разбиение осуществлялось. Соответственно, и множество модельных параметров с фиксированными значениями будет в общем случае разным в разных точках пространства модельных параметров. Это должно алгоритмически учитываться в процессе калибровки модели.
Рис. 1. Общая схема методологии численного прогнозирования характеристик природных динамических систем
6. Общая схема методологии численного решения Основной задачи
Рассмотрим теперь общую схему методологии, которая предлагается в данной работе в качестве базовой для численного решения Основной задачи. Схема представляет собой (рис. 1) совокупность выполняемых в определенной последовательности этапов. Рассмотрение общей схемы предлагаемой методологии проводится на уровне описания ее важнейших элементов, реализация которых является принципиально необходимой. Следует отметить, что сложная структурная организация реальных природных систем, многообразие и взаимное влияние протекающих в них процессов различного
6
Более подробно количественные оценки чувствительности модельных параметров обсуждаются в Приложении.
804
Вестник МГТУ, том 13, №4/1, 2010 г.
стр.793-809
характера предполагают необходимость использования развитого человеко-машинного взаимодействия как одного из важнейших элементов разрабатываемой методологии их численного анализа. Особенности организации такого взаимодействия обусловлены специфическими особенностями исследуемой природной системы и целевых задач проводимого исследования.
Обсудим основные этапы предлагаемой общей схемы методологии численного анализа природных систем с некоторыми необходимыми комментариями, безотносительно к конкретным природным системам и задачам. Будем предполагать при этом, что все требуемые действия выполнимы.
1. Сбор и обработка исходной информации о системе ДО и ее характеристиках. На этом этапе доступная априорная или экспертная информация о системе анализируется и подвергается необходимой предобработке с тем, чтобы появилась возможность задать все необходимые для исследования сведения. В том числе:
- временной отрезок [0,7];
- временную сетку п = {T1 / 7 = T1-1+ht1, hti>0, i=\,2,...,nt; Т0 = 0, Tnt = 7};
- необходимые сведения о внешних условиях формирования системы для каждого из временных интервалов 5T1=(T1-i,T1), i=l,2,...,nt;
- необходимые характеристики (стартовые) свойств слоев ft)i(0, i=1,2,...,nn(t) при вхождении их в исследуемую систему;
- конечное множество калибровочных разрезов Vk и дискретных полевых распределений P k(z) целевой характеристики для каждого калибровочного разреза Vk, £=1,2,...,nV при t=T - в настоящее время.
2. Построение £g(7) - структурной модели системы. Построение исходного варианта £g(7) -структурной модели системы с учетом динамики ее развития на временном отрезке [0,T]. Определение основных элементов структурной модели. Принятая при этом степень детализации их описания, осуществляются, как правило, экспертным путем на основании информации, заданной на первом этапе.
Следует отметить, что первый и второй этапы, будучи обязательными элементами методологии, не являются предметом исследования в данной работе (начальная структурная модель системы считается заданной). То же относится и к этапу 9.
3. Построение эффективной модели исследуемого процесса. Формирование эффективной математической модели исследуемого процесса и осуществление при необходимости ее декомпозиции. Определение перечня и задание диапазона возможных изменений значений модельных параметров. Этот этап полностью зависит от специфических особенностей конкретной системы, ее структурной модели и специфики исследуемого процесса, а также структуры полевых данных. Построение модели осуществляется в масштабе времени развития системы.
4. Постановка прямой задачи. Разработка численных методов ее решения. Постановка прямой задачи осуществляется в соответствии с выбранной эффективной математической моделью процесса с учетом структуры и особенностей полевых данных о значениях целевой характеристики систем. При построении численных методов решения прямой задачи необходимо учитывать ее математическую специфику, в том числе динамику развития исследуемой системы.
5. Постановка обратной задачи. Построение численных методов ее решения. Постановка обратной задачи не зависит от специфики исследуемой системы и может быть осуществлена в общем случае. Разработка метода решения обратной задачи осуществляется на основе регуляризирующего подхода, ориентированного на построение е-квазирешений.
Рассмотренные первые пять этапов общей схемы являются подготовительными этапами, предваряющими собственно процесс решения основной задачи. Все они выполняются при непосредственном участии специалистов в соответствующей предметной области. При этом третий, четвертый и пятый этапы выполняются с соблюдением изложенных выше основных принципов и, как правило, лишь один раз при первоначальной разработке методологии. При многократном практическом использовании методологии для численного анализа той или иной природной динамической системы эти этапы могут быть исключены из общей схемы.
6. Уменьшение количества модельных параметров. Упрощение используемых моделей посредством уменьшения количества модельных параметров. Осуществляется с помощью алгоритмических процедур, разрабатываемых в соответствии с рассмотренными выше идеями и ориентированных на анализ возможности:
- объединения соседних слоев (0,(1), ®1+1(t), i=1,2,...,nn(t)-1 в структурной модели системы Xg(0 в один слой;
- исключения из числа варьируемых модельных параметров тех параметров, фиксация значений которых не оказывает существенного влияния на качество калибровки модели исследуемого процесса.
805
Середа А.-В.И. Эффективный подход к численному моделированию...
Весь процесс организовывается, как правило, в интерактивном режиме. После завершения этого этапа формируется окончательный (рабочий) вариант структурной модели исследуемой системы и определяется рациональное количество варьируемых модельных параметров при решении обратной задачи.
7. Калибровка модели по имеющимся полевым данным. В зависимости от специфики структурной модели системы и пространственной структуры полевых данных, возможна различная организация выполнения данного этапа. Рассмотрим две характерных ситуации.
Первая из них имеет место в случае, когда по тем или иным причинам нецелесообразно проводить пространственную декомпозицию прямой задачи на одно- и двухмерные в пространственном отношении задачи, и калибровка модели осуществляется по всей пространственной области Q(T) в целом. В общем случае такая калибровка является очень трудоемкой в вычислительном отношении. Ситуация в практическом отношении нежелательная.
Вторая ситуация имеет место в случае, когда указанная пространственная декомпозиция прямой задачи целесообразна. Например, когда полевые данные пространственно структурированы и представлены лишь для конечного набора локальных вертикальных (калибровочных) разрезов Vkofl(T), £=1,2,.,nV. Тогда производится калибровка вертикальной математической модели для каждого калибровочного разреза Vk, £=1,2,. ,nV, которая осуществляется, например, посредством построения е-квазирешений обратной задачи (7), сформулированной в терминах задачи (18-21), с целью определения значений компонент соответствующих калибровочных векторов модельных параметров Х^*, £=1,2,...,nV. В зависимости от конкретных целей, калибровка может проводиться либо в автономном режиме для каждого калибровочного разреза, либо в режиме одновременной калибровки всех имеющихся калибровочных разрезов с дополнительными условиями на согласование получаемых результатов (региональная калибровка). Установка значений управляющих процессом констант осуществляется в интерактивном режиме.
8. Построение прогноза в заданной пространственной области. Этот этап выполняется полностью в автоматическом режиме. Построение прогноза осуществляется в два шага. На первом шаге посредством послойной интерполяции значений компонент калибровочных векторов Xk)*, £=1,2,.,nV в вертикальный разрез V0(x0,y0,z) определяются значения компонент вектора Х0)* - вектора модельных параметров для разреза V0. На втором шаге построение требуемого прогноза дискретного распределения P0(z) осуществляется в результате решения прямой задачи:
P0(z)= GZ(X0)*).
Здесь следует отметить, что послойная интерполяция значений компонент калибровочных векторов Xk)*, £=1,2,.,nV в вертикальный разрез V0(x0,y0,z) оправдывается сделанным ранее предположением о том, что пределах каждого слоя (01(t), i=1,2,...,na(t) структурной модели Xg(0 исследуемой системы для значений модельных параметров, как правило, характерна незначительная вариабельность. При этом в целом изменение этих значений в пространстве и во времени закономерно с физической точки зрения.
9. Анализ полученных результатов прогнозирования. Проводится экспертно. При неудовлетворительном заключении в интерактивном режиме может осуществляться дополнительное исследование имеющихся данных о системе и ее свойствах, целесообразная корректировка структурной модели системы Xg(t) и полевых данных, включая их возможное частичное исключение или дополнение, с целью повторения процесса прогнозирования, начиная с этапа 7.
10. Уточнение прогноза при поступлении дополнительной информации. Этот этап является дополнительным элементом предлагаемой методологии и может быть реализован при появлении различного рода дополнительной информации о свойствах исследуемой системы в окрестности прогнозного разреза V0. Уточнение прогноза целесообразно осуществлять при обнаружении значимого расхождения прогнозировавшихся и наблюдаемых значений прогнозируемой целевой характеристики.
7. Заключение
Рассмотренная общая схема методологии численного моделирования характеристик природных динамических систем в целом имеет общий характер. В методологическом плане она позволяет в полной мере реализовать сформулированную концепцию общего подхода к организации решения Основной задачи. Этот подход был успешно реализован применительно к актуальной практической проблеме численного прогнозирования значений геофлюидальных давлений7 в осадочных бассейнах земли (Мадатов, Середа, 2000а; 2000б; 2003; Madatov, Sereda, 2001; 2005а; 2005б; 2006).
7 давлений, воздействующих на флюид в осадочных толщах земли.
806
Вестник МГТУ, том 13, №4/1, 2010 г.
стр.793-809
Приложение. Чувствительность модельных параметров
Введем количественную оценку чувствительности модели в точке ХеХ0 к малому изменению значения модельного параметра /-1,2,.., т. Она может быть определена, например, как конечно-
разностное отношение:
S1(X,SXi)= (F(%+eSxi)"F(X,) )/$Х , /=1,2,..., т, (П1)
где F(x) - функция рассогласования, определенная в соответствии с (22); в1 - единичный вектор, все компоненты которого равны нулю, кроме компоненты с номером /, а размерность совпадает с размерностью вектора модельных параметров Xе X0; Sx1 - приращение значения компоненты вектора ХеХ0 Очевидно, что в этом случае величина S^Sx) эквивалентна значению конечно-разностной аппроксимации производной функции рассогласования F(x) по модельному параметру
Возможен и несколько другой подход к определению количественной оценки чувствительности модели в точке ХеХ0 к малому изменению значения модельного параметра X, /=1,2,..., т, когда указанная оценка вычисляется в соответствии с соотношением:
SXX)= l|G (X+e15Xi)-G(X)ii / l|G(X)ll, /=1,2,., m . (П2)
Очевидно, что использование формулы (П2) для определения значения S1(x,Sx1) обеспечивает большую вычислительную устойчивость по сравнению с использованием формулы (П1) при малых значениях Sx1 - приращения модельного параметра. Вместе с тем, в отличие от формулы (П1), формула (П2) всегда приводит к неотрицательной оценке чувствительности, что в определенных случаях может считаться недостатком.
Значение оценки S^Sx) будем называть чувствительностью модельного параметра X, /=1,2,.,т. Очевидно, что определенная указанным образом чувствительность модельного параметра X носит локальный характер и зависит, в общем случае, от калибровочного разреза Vq<^Q(T), q=1,2,...,nV, выбранной точки XeX0 и величины Sx.
Для того, чтобы в определенной степени повысить объективность оценки как для отдельного разреза, так и на всем множестве калибровочных разрезов (в области £2(T) в целом), а также сопоставимость таких оценок для различных модельных параметров и разрезов, предлагается использовать следующие правила при их вычислении.
Во-первых, будем использовать так называемое правило нормировки. А именно:
- диапазон изменения всех модельных параметров должен быть приведен к отрезку [0,1], что осуществляется в результате простого масштабирования в соответствии с правилом:
X1:=(X1-X1 ш1п)/(Х1 max -X mmX ^У^.^ m; (П3)
- величину 8Xi - приращение модельного параметра X, используемую для расчета оценки чувствительности, примем одинаковой для всех /=1,2,.., т.
X = X /=1,2,.., т. (П4)
Например, Sx = 0.01. В дальнейшем Sx будем называть также нормированным приращением модельных параметров.
При проведении практических расчетов значение Sx задается экспертным путем или определяется в результате вычислительного эксперимента. При этом главным критерием выбора значения Sx является его максимально возможное значение, при котором имеет место определенная стабильность получаемых значений S1(x,Sx1) - оценок чувствительности для всех модельных параметров задачи. Основанием для такого выбора является постулируемое свойство локальной однородности пространства модельных параметров в достаточно малой окрестности заданной точки. Подобное предположение является вполне естественным для характеристики свойств физической среды, в которой протекают изучаемые процессы. В данном случае оно распространяется на используемую математическую модель этих процессов и заключается в предположении, что свойства математической модели считаются практически неизменными в достаточно малой окрестности любой точки пространства модельных параметров. В том числе и такое свойство модели, как ее чувствительность.
Во-вторых, будем следовать принципу или правилу определения оценки чувствительности модели в соответствии с результатом ее калибровки (калибровочному правилу), согласно которому расчет чувствительности модельных параметров осуществляется не в произвольной точке хеХ°, а в точке Ха еХ0, соответствующей решению обратной задачи для конкретного калибровочного разреза. Такой подход вполне оправдан, поскольку калибровочные векторы Ха, полученные из условия согласования
807
Середа А.-В.И. Эффективный подход к численному моделированию...
полевых и модельных данных, в эффективном смысле отражают свойства исследуемой системы в окрестности калибровочных разрезов, и анализ чувствительности модели в этих точках представляется наиболее адекватным в отношении изучаемых процессов.
Основываясь на сформулированных выше правиле нормировки и калибровочном правиле для каждой скважины VqczQ(T), q=1,2,...,nV по формулам (П1) или (П2) может быть осуществлен расчет чувствительности по всем модельным параметрам %1, /=1,2,.., т.
В результате для каждого калибровочного разреза Vq будет сформирован Sq(%qE ,8%) - вектор чувствительности модели, компонентами которого являются чувствительности соответствующих модельных параметров. Компоненты вектора чувствительности Sq(%qE,8%) модели могут быть наглядно представлены в виде диаграммы чувствительности. При ее построении будем следовать для удобства двум условиям. Во-первых, все оценки чувствительностей модельных параметров неотрицательны. Если это на самом деле не так, что может быть при использовании для их вычисления правила (П1), они берутся по модулю. Во-вторых, все значения оценок чувствительностей нормируются таким образом, чтобы они принадлежали отрезку [0,1]. Пример такой диаграммы приведен на рис. П1. Линии определенного цвета отображают на этой диаграмме чувствительности соответствующего формационного модельного параметра.
Анализ компонент векторов Sq(%qE ,8%) как для одного отдельного разреза, так и для всех разрезов VqcQ(T), q=1,2,...,nV в совокупности является одним из важных элементов рассматриваемых ниже подходов к построению эффективных моделей с меньшим числом модельных параметров8.
13 -
12 . ...............
.И ....... ........... ~ “
Рис. П1. Диаграмма чувствительности
Нормированная чувствительность модельных параметров
Классификация калибровочных разрезов
Очевидно, что получаемые в результате калибровки вертикальной задачи калибровочные векторы %qE, q=1,2,...,nV не являются единственным решением соответствующей обратной задачи. Существует бесконечное множество векторов %qeX°, являющихся ее 5е-квазирешениями, то есть удовлетворяющих условию:
F(Xq) < 4=8е2, %qeX°. (5)
Эти решения неразличимы с точки зрения (П5) и являются в этом смысле равноценными решениями обратной задачи. Они принадлежат некоторому подмножеству XqE0 множества X0, для всех точек которого выполнено условие (П5). В общем случае при достачочно малом значении (а также при неудачном задании множества X°) XqE0 может представлять собой объединение s не пустых попарно непересекающихся односвязных подмножеств XqEJ°, j=1,2,.,s:
XqE0 = Xq/U XqE2°U X^U . . . U Xq^, (П6)
где
XqEj ^0,j=1,2,...,s; XqEJnXqE1=0, V tej, /'=1,2,...s; j=1,2,...s. (П7)
8 Анализ чувствительности модельных параметров оказывается полезным и в процессе калибровки модели.
808
Вестник МГТУ, том 13, №4/1, 2010 г.
стр.793-809
В дальнейшем будем предполагать, что для каждой калибровочной скважины Vq^Q(T), q=1,2,...,nV значение 5=1. Это означает, в частности, что для всех калибровочных скважин XqE0rf° является односвязной (в силу непрерывности Gq(x) - оператора прямой задачи) подобластью X0. Данное предположение заметно упрощает последующие рассуждения, не умаляя при этом их общности (уменьшая длины интервалов возможного изменения модельных параметров всегда можно соответствующим образом задать множество Х°). В этом случае будем считать, что получаемое в результате решения обратной задачи ее д£ -квазирешение xqe единственным.
В этом контексте возможно осуществление классификации (группировка) калибровочных разрезов VqcQ(T), q=1,2,...,nV по признаку возможной принадлежности соответствующих им калибровочных векторов %qs одному и тому же подмножеству XE°cX°, где XE° - односвязная подобласть X0, представляющая собой непустое пересечение всех подмножествXqs0cX°, q=1,2,...,nV:
Xs° = Xi E°n X2E°n XsE°n ... n XnvEV0. (8)
Признаки возможности существования подмножества XE°, удовлетворяющего (П8),
устанавливаются, как правило, экспертным путем. К формальным признакам такого рода можно отнести, например, достаточную близость значений одноименных компонент калибровочных векторов для разрезов одной группы, либо их физически интерпретируемую изменчивость.
В качестве заключительного замечания отметим, что пространственная "конфигурация" подмножества XE°cX° совершенно определенным образом зависит от чувствительности модели к изменению значений различных модельных параметров. Подмножество XE° "растягивается" вдоль осей, соответствующих нечувствительным параметрам, и сжимается вдоль осей, соответствующих чувствительным модельным параметрам. Указанная пространственная "конфигурация" XE° может "исправляться" посредством уменьшения длины интервалов возможных значений нечувствительных параметров вплоть до "фиксирования" этих значений.
Литература
Hadamard I. Le probleme de Cauchy et les equations aus derivees partielles lineaires hyperboliques. Hermann, 1932.
Madatov A.G., Sereda V.-A.I. A multi-well data inversion purpose-built for calibration of an effective basin model at overpressure prediction. Proceedings of the Murmansk State Technical University, v.8, N 1, p.44-83, 2005б.
Madatov A.G., Sereda V.-A.I. The decomposition of 3-D overpressure evolution model in basin scale and its application to the fault seal analysis. Proceedings of the Murmansk State Technical University, v.4, N 1, p.79-96, 2001.
Madatov A.G., Sereda V.-A.I. The effective basin model concept and fast 3-D overpressure modeling in basin time scale. Proceedings of the Murmansk State Technical University, v.8, N 1, p.5-43, 2005а.
Madatov A.G., Sereda V.-A.I. The pre drill and real time while drilling overpressure prediction based on Effective Bain Model concept. Proceedings of the Murmansk State Technical University, v.9, N 3, p.361388, 2006.
Алифанов О.М., Артюхин Е.А., Румянцев С.В. Экстремальные методы решения некорректных задач. М., Наука, 288 с., 1988.
Бакушинский А.Б., Гончарский А.В. Итеративные методы решения некорректных задач. М., Наука, 128 с., 1989.
Гончарский А.В. Некорректно поставленные задачи и методы их решения. В сб.: Некорректные задачи естествознания. Под ред. А.Н. Тихонова, Л.В. Гончарского. М., Изд-во МГУ, c.15-36, 1987.
Иванов В.К. О линейных некорректных задачах. ДАН СССР, т.196, № 145, с.270-272, 1970.
Иванов В.К. О некорректно поставленных задачах. Математ. сб., т.61, № 2, c.211-223, 1963.
Мадатов А.Г., Середа В.-А.И. Прямая и обратная задача геофлюидодинамики в приложении к прогнозированию зон АВПД в осадочных бассейнах. 1. Теоретический аспект. Вестник МГТУ, т.3, № 1, c.89-114, 2000а.
Мадатов А.Г., Середа В.-А.И. Прямая и обратная задача геофлюидодинамики в приложении к прогнозированию зон АВПД в осадочных бассейнах. 2. Практический аспект. Вестник МГТУ, т.3, № 2, c.351-366, 2000б.
Мадатов А.Г., Середа В.-А.И. Рационализация уровня сложности бассейновой модели среды для целей прогнозирования свойств геофлюидальной системы. Вестник МГТУ, т.6, № 1, c. 119-144, 2003. Тихонов А.Н. Об устойчивости обратных задач. ДАН СССР, т.39, № 5, c.195-198, 1943.
809
