Вариационное усвоение данных в реальном времени Текст научной статьи по специальности «Математика»

ВАРИАЦИОННОЕ УСВОЕНИЕ ДАННЫХ В РЕАЛЬНОМ ВРЕМЕНИ*

В. В. ПЕНЕНКО Институт вычислительной математики и математической геофизики СО РАН, Новосибирск, Россия

e-mail: [email protected]

A brief description of the lecture notes on data assimilation given to the participants of the Young Scientists School on Computational Information Technologies for Environmental Sciences "CITES-2005" (Novosibirsk, March 13-23, 2005) is provided. New original algorithms of successive data assimilation are presented. The decomposition and splitting techniques used in these algorithms are based on variational principles. One such algorithm was proposed to the students for practical training.

Введение

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

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

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

* Работа выполняется по Программам фундаментальных исследований РАН (грант № 13), ОМН РАН (грант 1.3.2), поддержана Российским фондом фундаментальных исследований (грант № 04-05-64562) и Интеграционными грантами СО РАН (№ 03-130, № 03-137, № 03-138).

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2005.

В работе [9] дан пример применения этой методики для реализации 4-мерного усвоения данных на базе глобальной модели динамики атмосферы в рамках метода расщепления. Как отмечено в обзоре [10], это были первые в мире работы по вариационному усвоению данных в метеорологии с использованием сопряженных уравнений. Большой вклад в теорию и практику вариационного усвоения внесла французская школа [10, 11]. Новый этап в разработке вариационных алгоритмов начался с идеи формулировки задачи усвоения с учетом неопределенностей (ошибок) в моделях и данных наблюдений [12, 13]. В этом случае вариационный принцип строится для минимизации функционалов, содержащих суммарную меру всех неопределенностей. Интересно отметить, что алгоритмы типа Кал-мана можно получить из алгоритмов с сопряженными уравнениями. В частности, было показано, что в случае линейных моделей с квадратичными функционалами качества схемы с сопряженными уравнениями и алгоритмы типа фильтрации Калмана эквивалентны, т. е. представляют собой две различные схемы реализации условий стационарности одних и тех же функционалов. В нелинейных случаях они эквивалентны по точности по отношению к целевому функционалу.

Поскольку усвоение данных уже вошло в оперативную практику, особое внимание в настоящее время уделяется созданию таких вариационных алгоритмов, которые могли бы работать в реальном времени [14, 15, 19-21]. Много работ в последнее время посвящено усвоению химических и аэрозольных данных [16], имеющих непосредственное отношение к качеству окружающей среды.

В настоящей статье кратко изложены основные положения курса лекций по усвоению данных, прочитанного автором слушателям Международной школы молодых ученых по вычислительно-информационным технологиям для наук об окружающей среде "С1ТЕЯ-2005" (Новосибирск, 13-23 марта, 2005). Представлены новые оригинальные алгоритмы последовательного усвоения данных, использующие методы декомпозиции и расщепления в рамках вариационных принципов. Один из таких алгоритмов был предложен участникам школы для реализации в качестве компьютерной практики.

1. Модели процессов

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

Ц<p) = B^ + G(<p,Y)-f-т = 0, (1)

где В — диагональная матрица; V) — "пространственный" нелинейный матричный

дифференциальный оператор, в числе основных элементов в нем — адвективно-диффузи-онные операторы, действующие на различные компоненты функций состояния; ¥ — функции источников; г — функции, описывающие неопределенности моделей. При необходимости операторы переноса симметризуются с помощью уравнения неразрывности [3].

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

Здесь = {</?Дх, ¿),г = 1, гг} € (^(И^ — вектор-функция состояния, ее компоненты описывают потенциальную температуру, отношения смеси для характеристик влажности в атмосфере (водяной пар, облачная вода, дождевая вода, снег и ледовые кристаллы), концентрации загрязняющих примесей в газовом и аэрозольном состоянии; ^ = {/(х, £), г =

1, п } — функции источников тепла, влаги и примесей; гг — функции, описывающие неопределенности и ошибки моделей; п — функция, зависящая от выбора системы координат; и = (и1,и2,и3) — вектор скорости; = diag{^1, — коэффициенты турбулентного обмена для субстанции в направлении координат х = {ж», г = 1, 3} £ 1), Д = Д х [0, (Н(ф))г — нелинейные матричные операторы, которые описывают локальные процессы трансформации соответствующих субстанций. Все необходимые элементы гидродинамического фона для моделей (2) рассчитываются по соответствующим моделям динамики атмосферы, согласованным с (2) в рамках вариационного принципа на уровне интегральных тождеств и их сумматорных аналогов. Функции п, и, ¡лг, /г, внутренние параметры операторов Нг, входные данные начальных и краевых условий для моделей (2) включаются в совокупность компонент вектора параметров V.

2. Вариационная формулировка моделей процессов

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

Структура областей и систем координат, основные обозначения, постановки задач, краевых и начальных условий, описание функций состояния и параметров, структура функционалов в зависимости от целей исследований и способы построения тождеств для вариационных формулировок моделей описаны, например, в [3, 4].

Наряду с дифференциальной постановкой задачи будем использовать вариационную формулировку модели (1), (2)

где £ Q* — функции, сопряженные по отношению к Интегральное тождество (3) строится с учетом краевых и начальных условий так, чтобы при подстановке ^* = ^ функционал (3) давал соотношение баланса энергии исследуемой системы, а при = const —

(3)

соотношения баланса первого порядка. Например, выполнив все необходимые преобразования в (3) для модели (2), окончательно получим интегральное тождество в виде

IУ, <р*) = ¿| (Л^, Л + /((Я- ¡г - т№пйВйЛ = 0. (4)

Здесь

(Л^*)г = 0.5

</? --"гулч/Ш — </зшу7г(/3 и)

+

grad^ grad^*j>dDdt ^ J 0.5)^^*пdD |ц +

Б

+ J (о.Ь(рип-I((В.ьч>) - q)ч>*'кd£ldt} , (5)

Пг Пг

ип — нормальная к границе составляющая вектора скорости. Заметим, что операторы Яь в граничных условиях для системы (2) желательно записать в таком виде, чтобы можно было с их помощью исключить слагаемые с нормальной производной от функции состояния в двух последних интегралах в (5).

Формы (Л^, представляют операторы переноса и турбулентного обмена. В них операторам переноса соответствуют антисимметричные относительно функций ^ и пары слагаемых, а операторам турбулентности — симметричные слагаемые. Свойства антисимметрии и симметрии в принципе должны сохраняться при аппроксимации интегралов и подынтегральных выражений. Кстати, отметим, что замену переменных и систем координат проще проводить в рамках интегрального тождества.

Для согласованного расчета динамики атмосферы и переносов примесей в тождество аддитивно включаются слагаемые, соответствующие уравнениям движения [3].

3. Модели наблюдений и функционалы

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

Фт + ^(М), (6)

где Фт — набор наблюдаемых величин; W(^) — совокупность моделей наблюдений; п(х, £) — ошибки и неопределенности этих моделей и данных. Значения Фт определяются на множестве точек D1jn € Dt. Символ [ ]т обозначает операцию переноса информации с Dt на Dtm, например, с помощью операторов проектирования или интерполяции.

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

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

где индекс Т означает транспонирование; Б — весовая матрица для формирования скалярного произведения на множестве данных наблюдений; Хо — весовая функция, определяющая конфигурацию пространственно-временного носителя наблюдений в ^ и меру для представления (7) в виде соответствующих интегралов по области Dt. Для обнаружения источников и планирования наблюдений, в дополнение к функционалам (7), необходимо рассмотреть последовательность функционалов, каждый из которых описывает индивидуальные наблюдения в составе (6).

Традиционный подход к решению задач охраны окружающей среды обычно базируется на методах прямого моделирования. Суть их состоит в создании математических моделей изучаемых процессов и в проведении сценарных расчетов при различных способах задания входных данных и внешних воздействий. Несмотря на широкое распространение этого подхода, он не может обеспечить решение всего комплекса возникающих вопросов на современном уровне. Многие задачи требуют комбинированного использования методов прямого и обратного моделирования [3, 8, 12, 13, 19]. В таком варианте прямое моделирование представляет модели процессов различной степени сложности по функциональному содержанию и детализации, по степени учета различных факторов и пространственно-временному разрешению дискретных аппроксимаций. Методы обратного моделирования реализуют более высокие системные уровни проблемы в целом, в которых модели процессов выступают в качестве связей между функциями состояния, входными параметрами и источниками внешних воздействий. Организация такой технологии моделирования строится на базе вариационных принципов в сочетании с методами декомпозиции, расщепления и комплексирования [3, 8, 13]. Перечислим некоторые типовые задачи, решаемые с помощью комбинации методов прямого и обратного моделирования: диагностика качества моделей; усвоение данных измерений; комплексирование моделей различных масштабов; восстановление функций состояния; исследование чувствительности моделей к вариациям входных данных; планирование наблюдений и оценки информативности систем мониторинга; оценки наблюдаемости территорий; экологическое проектирование с позиций устойчивого развития; оценка мощности и локализация местоположения источников; управление источниками; районирование территории в соответствии с уровнями антропогенных нагрузок; оценка областей влияния и опасности источников; задачи типа "рецептор — источник — рецептор"; оценки риска и уязвимости по оношению к антропогенным воздействиям.

4. Основной алгоритм обратного моделирования

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

(7)

проксимации получаются с учетом структуры функционала качества и способов нахождения его стационарных значений на множествах значений функций состояния, параметров моделей в дискретной формулировке [3, 8, 13].

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

Ф £ (V) = +

т> N £

0\Т ^ (,„0

(пТМ,П';,„ + (гТМ2г)^ + ((V0 - чИ? Мз (V - +

К

+ ((У - У«)Т М4 (У - Уа)) ^ л /2 + У, V*), к > 1. (8)

Здесь первое слагаемое представляет собой целевой функционал, следующие четыре члена суммарно выражают меру неопределенностей моделей наблюдений, моделей процессов, начальных данных и параметров соответственно. Они строятся по аналогии с (7). Последнее слагаемое содержит описание численной модели в вариационной формулировке (3)-(5); Мг, г = 1,4, — весовые матрицы; Уа и — априорные оценки.

Дискретизацию функционалов (8) осуществляем так же, как и (3)-(5), с использованием методов слабой аппроксимации, расщепления и декомпозиции. Все кубатурные формулы для интегралов строятся на одних и тех же сетках. Окончательно дискретные аппроксимации моделей и алгоритмы моделирования получаются из условий стационарности функционалов Ф

1) для основных задач и методов прямого моделирования — из условий стационарности к вариациям компонентов сопряженной функции V*;

2) для сопряженных задач и методов обратного моделирования — из условий стационарности к вариациям компонентов функции состояния V;

3) если в моделях учитываются функции неопределенностей, то условия стационарности к вариациям компонентов этих функций дают систему уравнений для их расчета по фактической информации, заложенной в функционалах в составе (8) через соответствующие функции чувствительности;

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

В результате получается следующая система операторных уравнений [13]:

дФ£

^ = ВА& + Gh((p,Y)-f-т = 0] (9)

дФ £

= (ВА<)т<р*к + Ат(<р, У)<р*к + а, = 0; (10)

^к(х) = 0; (11)

д

= —(ф^ + О^М^)); (12)

V0 = V + Мз-1 VI (0), г = 0; (13)

г(х,г) = М2-^к(х, г); (14)

У = Уа + М-1 Гк; (15)

д

= (16)

д г п

У)<р' = — [С\ц> + а<р', У)] |а=0 . (17)

Здесь Л^ — оператор производных по времени или их дискретных аппроксимаций; Ат(ф, У) — пространственный оператор сопряженной задачи; Г& — функции чувствительности моделей к вариациям параметров; а — вещественный параметр; р = — вариации функции состояния. Операции дифференцирования в (9), (10), (12), (16), (17) осуществляются для всех сеточных компонентов функции состояния, сопряженной функции и параметров. Конструктивно они реализуются с помощью производных Гато для функционалов (4), (5), (8) относительно всех их функциональных аргументов в дискретном представлении. Выбранный способ дискретизации функционалов приводит к тому, что полученные с помощью вариационного принципа дискретные аналоги основных (9) и сопряженных (10)-(12) задач представляют собой взаимно согласованные схемы расщепления и системы моделей в рамках выбранного способа декомпозиции.

В сопряженных задачах градиенты функционалов (12) относительно компонентов функции состояния в узлах сеточной области выступают в качестве функций источников при организации процедур усвоения данных дистанционных и контактных наблюдений и при учете с помощью целевых функционалов ограничений в оптимизационных задачах управления и проектирования. Так как функции состояния входят практически во все слагаемые в (8), сопряженные задачи по структуре и функциональному содержанию замыкают на себя все внутренние связи между различными элементами системы моделирования, учитываемыми в основном функционале. При этих условиях оценки вариаций целевых функционалов получаются оптимальными в том смысле, что они имеют точность второго порядка малости относительно вариаций ф, ф*, г. Система уравнений (9)-(17) решается относительно ф, ф*, г, ф0, У итерационными процедурами, начиная с г = 0, ф0 = фа, У = Уа. Эта система представляет собой центральное ядро вычислительной технологии прямого и обратного моделирования для решения перечисленных выше задач и их модификаций.

5. Алгоритмы фильтрации Калмана

Если в алгоритме (9)—(17) исключить функцию ф*, то получаются схемы усвоения типа фильтрации Калмана:

^ + У)-?- РитМЛЧ - Шу?)) = 0; (18)

т

^ц=о = V"0

дР

фк=о = фа; (19)

+ АР + РАТ + ритм1ир - М^1 = 0; (20)

и^Щ^, Р(х,0) = ^(х>. (21)

Уравнение (18) представляет собой модификацию уравнения (9). Здесь Р = Р(х,{) — весовая матрица фильтра Калмана. Ее размерность равна квадрату размерности векторов

ф в дискретном представлении на сетке Она вводится для включения невязки п из (6) непосредственно в уравнение прямой модели (9) вместо функции г, которая учитывает эти невязки посредством функций ф* с помощью (13),(14). Уравнения (20),(21) для нахождения матриц Р как раз и получаются из требований исключения функции ф* из

(9),(13),(14) .

За почти 50-летнюю историю развития методов фильтрации Калмана разработано множество их модификаций. На практике часто используются дискретные варианты адаптивного усвоения (см., например [17, 18]). Следует подчеркнуть, что хотя эти алгоритмы выведены из других соображений, их также можно построить исходя из вариационного принципа, вытекающего из метода взвешенных наименьших квадратов (МВНК). Алгоритмы этого типа реализуются в два этапа.

1. На каждом временном интервале < £ < решается прямая задача (9) с г = 0 при условиях (рЭ-1 = Ср3~ , ^ = 1, Jkl ф = ф°а- В результате получается решение фК

2. Далее осуществляется корректировка решения с учетом данных наблюдений. Для этого определяется целевой функционал МВНК с использованием обозначений (8)

(Ф(Ф)У = \{(<Р- ф)ТМ3(ф - ф) + (Ф - Ш(ф))тМ1(^ - \¥(<р))У , (22)

где ф — искомая функция состояния; М3, М1 — заданные положительно определенные весовые матрицы. Индекс ] указывает, что все объекты рассматриваются на шаге ,tj]. Искомое решение = фj получается из условия стационарности (22) к вариациям

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

ф = {( + К(Ф - W(())} ; (23)

Kj = {РитМ1}5 , Р5 = | [М3 + итМ1и1. (24)

Матрицы К и Р, вычисляются при ф = ф. Если данные наблюдений отсутствуют, то полагается ф5 = ф.

Для нелинейных моделей процессов и моделей наблюдений схема (18)-(21) имеет точность 0($(2), где 5( = Р(*, по отношению к схеме (9)-(17). А это значит, что обе схемы эквивалентны по точности. В случае линейных операторов моделей обе схемы полностью эквивалентны и система (18)-(21) является классической схемой типа Калмана [6].

Таким образом, алгоритмы с сопряженными уравнениями и алгоритмы типа фильтрации Калмана имеют одну и ту же идейную основу, выраженную вариационным принципом для оценок функционала типа (8). Конструктивно они реализуют разные пути для отыскания стационарных состояний целевых функционалов. Сопряженные функции в алгоритмах (9)-(17) замыкают на себя все внутренние степени свободы в соответствии с целевыми функционалами (8). В отличие от этого, в алгоритмах типа фильтрации Калмана роль замыкающего звена играют весовые матрицы Р и К. При большой размерности задач (функции состояния в современных численных моделях имеют размерность 106 ... 108) табулирование этих матриц для рассматриваемого класса моделей представляет собой трудновыполнимую сверхзадачу как по объемам вычислений, так и с позиций устойчивости. Этим обусловлено появление большого числа облегченных модификаций, так называемых "субоптимальных" алгоритмов, оценить точность которых по отношению к исходной задаче в общем случае не представляется возможным [18]. На практике применимость таких алгоритмов проверяется экспериментально.

6. Аддитивные алгоритмы усвоения данных

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

зк-1 зк-1

= £ ^; ^ = х [^,13], Ф?(ф, ф*, У, Ф) = £ Ф?(ф, ф*, У, Ф), (25) ¿=1 ¿=1

где Ф? — часть общего фукционала (8), относящаяся к интервалу времени ], ] =

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

Оставаясь в рамках общей схемы расщепления, реализованной в системе (9)-(17), рассмотрим ее модификации для усвоения данных в реальном времени. Для этого предположим, что в схеме расщепления имеются 5 этапов (в > 1). Первые 5 — 1 этапов в формуле (9) реализуются в обычном режиме прямого моделирования, а последний этап модифицируем исходя из целей усвоения данных. Запишем его в операторном виде

Л„ ф — % — 4 = 0, (26)

где Лsjф — оператор, аппроксимирующий часть модели в-го этапа расщепления; ^ — функции источников; г^ — функция, описывающая неопределенности модели в исходной постановке (1), (2) и дополнительные неопределенности, вносимые в численную модель схемой расщепления на шаге декомпозиции [tj-1,tj].

Рассмотрим случай, когда производятся измерения компонентов функции состояния, а точки наблюдения совпадают с некоторыми узлами сеточной области Б™ € Б?. Весовые матрицы M1j и M2¿ для оценок неопределенностей моделей и наблюдений в функционалах (8), (25) будем задавать диагональными. При таких предположениях определим локально сопряженную по отношению к (26) задачу из состава (10) и соотношения для оценок неопределенностей моделей (14)

Л^. ф^' = ау M1¿• (Ф^1 — ф^1); (27)

г2 = (М^/а^ )ф*. (28)

Здесь для упрощения схемы оценка невязки между измерениями и результатом моделирования берется по информации в момент времени t = tj-1, которая является исходной для

прогноза на интервале [tj-1,tj]. Полученная схема реализуется в следующем порядке: решаются задачи (27), (28), (26) по всем элементам схемы расщепления на рассматриваемом этапе. Заметим, что этот алгоритм был предложен участникам школы для реализации в качестве компьютерной практики по методам усвоения данных.

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

Л^. ф= ау Му (Фу — фу); (29)

Лу фу — / — (М- /ау )ф*у = 0. (30)

В этой схеме усваиваемая информация поступает формально в момент времени t = tj. Задача (28)-(30) представляет собой систему дискретных уравнений, которая решается прямым алгоритмом. Для моделей переноса (2) конечно-разностный оператор Л5 обычно представляется в виде трехдиагональной матрицы. В этом случае задачи (26), (27) решаются с помощью прямого метода трехточечной прогонки, а система (29), (30) — трехточечной матричной прогонки с матрицами второго порядка. Устойчивость вычислений по схемам (26)-(30) для задач переноса обеспечивается свойством диагонального преобладания в матрицах Лу и Лу. В общем случае устойчивость этих двух схем усвоения будет обеспечена, если схемы расщепления в (9), (10) устойчивы в целом. Такая устойчивость гарантируется свойством энергетической сбалансированности, заложенной в тождествах (3)-(5) и в их аппроксимациях в функционалах (8).

Приведем еще один тип схем усвоения данных в реальном времени [19, 20], которые получаются непосредственно в результате минимизации частей функционала качества, относящихся к шагу декомпозиции. Схемы расщепления типа (19), (26) используются в качестве ограничения на искомые функции состояния. Запишем функционал качества для этой схемы в виде

Фоу(ф) = 0.5 [а1 (пТщ) + а2(гТWljг у)] ; (31)

Г = ЛУ <Ру — Л, Пу = Фу — <Ру. (32)

Здесь а1, а2 > 0 — параметры, удовлетворяющие условию а1 + а2 = 1. Подставляя (32) в (31) и записывая условия минимизации полученного функционала относительно функции ф, запишем уравнение

л.:;л/|; - Г) + ^-М^ - Ф') = 0. (33)

а2у

Легко показать, что система (33) однозначно разрешима. Для уравнений переноса с трехточечной аппроксимацией адвективно-диффузионных операторов и с диагональными положительно определенными весовыми матрицами в функционале (31) матрица системы (33) получается пятидиагональной. Полученная система эффективно решается методом пятиточечной прогонки. Непосредственными вычислениями получаем, что схема (29), (30) эквивалентна (33), т.е. это два варианта реализации условий минимизации одного и того же функционала (31) при условии (26). Если в (33) невязку во втором слагаемом рассчитывать по априорным оценкам результатов наблюдений и функции состояния в момент времени t = t¿-1, то из (33) получается схема усвоения, эквивалентная (26)-(28).

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

расщепления и декомпозиции они имеют высокоэффективные схемы реализации для рассматриваемого класса задач. В этом плане им нет равных среди других схем усвоения данных. По точности все рассмотренные здесь схемы эквивалентны, поскольку они порождаются условием стационарности одного и того же функционала качества (8) в дискретном представлении. Функции неопределенностей, оцениваемые в алгоритмах по формулам (14), (28), дают дополнительную информацию для планирования наблюдений. Естественно предположить, что дополнительные наблюдения следует производить там, где функции неопределенности имеют относительно большие значения.

В случае более сложных моделей наблюдений (включая нелинейные) необходимо соответственно сформировать правые части в схемах (27), (29) и во втором слагаемом (33) по формулам типа (12), выражающим условия минимизации функционала неопределенности наблюдений в (8) и в его декомпозированных частях в (25) по отношению к компонентам функции состояния.

Список литературы

[1] Марчук Г.И. Математическое моделирование в проблеме окружающей среды. М.: Наука, 1982.

[2] Марчук Г.И. Сопряженные уравнения и анализ сложных систем. М.: Наука, 1992.

[3] Пененко В.В. Методы численного моделирования атмосферных процессов. Л.: Гидроме-теоиздат, 1981.

[4] Пененко В.В., Алоян А.Е. Модели и методы для задач охраны окружающей среды. Новосибирск: Наука, 1985.

[5] SASAKI I. An objective analysis based on variational method //J. Met. Soc. Japan. 1958. Vol. 36, N 3. P. 29-30.

[6] KALMAN R.E. A new approach to linear filtering and prediction problems // Trans. AME. J. Basic Eng. 1960. Vol. 82. P. 34-35.

[7] KALMAN R.E., BuOY R.S. New results in linear filtering and prediction theory // Trans. AME. Ser. D. J. Basic Eng. 1961. Vol. 83. P. 95-107.

[8] Пененко В.В. Вычислительные аспекты моделирования динамики атмосферных процессов и оценки влияния различных факторов на динамику атмосферы // Некоторые проблемы вычислительной и прикладной математики. Новосибирск: Наука, 1975. С. 61-76.

[9] Пененко В.В., Образцов Н.Н. Вариационный метод согласования полей метеорологических элементов // Метеорология и гидрология. 1976. № 11. С. 3-16.

[10] TALAGRAND O., Courtier P. Variational assimilation of meteorological observations with the adjoint vorticity equation. I: Theory // Quart. J. Roy. Meteor. Soc. 1987. Vol. 113. P. 1311-1328.

[11] Le Dimet F., TALAGRAND O. Variational algorithms for analysis and assimilation: theoretical aspects // Tellus. 1986. Vol. 38A. P. 97-110.

[12] Пененко В.В. Системная организация математических моделей для задач физики атмосферы, океана и охраны окружающей среды. Новосибирск, 1985. (Препр. ВЦ СО РАН; № 619). 43 с.

[13] Penenko V.V. Some aspects of mathematical modeling using the models together with observational data // Bulletin of the NCC. Series: Num. Model. in Atmosph. 1996. N 4. P. 32-51.

[14] DAESOU D., Carmiohael G. An adjoint sensitivity method for the adaptive location of the observations in air quality modeling //J. Atmos. Sci. 2003. Vol. 60, N 1. P. 434-449.

[15] TREVISAN A., Uboldi F. Assimilation of standard and targeted observations within the unstable subspace of the observation-analysis-forecast sycle system //J. Atmos. Sci. 2004. Vol. 65, N 1. P. 103-113.

[16] Elbern H., Sohmidt H., Talagrand O., Ebel A. 4D-variational data assimilation with an adjoint air quality model for emission analysis // Environmental Modelling and Software. 2000. Vol. 15(6-7). P. 539-548.

[17] БраЙСОн А., Хо Ю-Ши. Прикладная теория оптимального управления. М.: Мир, 1972. 544 с.

[18] Сейдж Э.П., Уайт Ч.С. Оптимальное управление системами. М.: Радио и связь, 1982. 392 с.

[19] Tsvetova E.A., Penenko V.V. Adjoint equations and splitting technique for fast data assimilation and sensitivity studies // The use of Data Assimilation in Coupled Hydrodynamic, Ecological and Bio-geo-chemical Models of the Ocean. Abstracts of the 33rd Intern. Liege Colloquium on Ocean Dynamics. 2001. P. 50-51.

[20] Penenko V.V., Tsvetova E.A. Variational Fast Data Assimilation Algorithms. Research Activities in Atmospheric and Oceanic Modeling. WGNE Blue Book Web, 2002. Site http://www. cmc.ec.gc.ca/rpn/wgne 01-48

[21] Пененко В.В., Цветова Е.А. Быстрое усвоение данных в атмосферных и океанических исследованиях // Вычисл. технологии. 2002. Т. 7. Специальный выпуск. С. 141-146.

Поступила в редакцию 21 июня 2005 г.