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

Вычислительные технологии

Том 11, часть 2, Специальный выпуск, 2006

НЕКОТОРЫЕ ТЕОРЕТИЧЕСКИЕ И ПРИКЛАДНЫЕ ВОПРОСЫ ПОСЛЕДОВАТЕЛЬНОГО ВАРИАЦИОННОГО

УСВОЕНИЯ ДАННЫХ*

А. В. Пененко Институт вычислительной математики и математической геофизики СО РАН, Новосибирск, Россия e-mail: aleks@ommgp.sscc.ru

A method for sequential variational data assimilation is considered. The Kalman filter is proved to be a special case of the algorithm. The economical algorithm for the state estimation of the 2D distributed system when the state function is measured at particular points is constructed.

Введение

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

1. Постановка задачи

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

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 04-05-64562), Программы Президиума РАН (№ 16), Программы Отделения математических наук РАН (№ 3), Европейской комиссии (контракт № 013427).

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

состояния щк € М", Также имеется некоторая гипотетическая модель системы с параметрами Лк+1 (вещественной невырожденной матрицей п х и), /к+1 € М", которая позволяет спрогнозировать па основе текущего состояния модели ее следующее состояние Щк+1 по алгоритму

Лк+1Щк+1 = щк + т/к+1. (1)

Пусть также есть данные некоторых измерений фк+1 € Мт, Если Щк+1 — реальное состояние системы, то пусть связь между ними задается как

Н к+1Щк+1 = фк+1 + пк+1,

где пк+1 — неизвестные ошибки измерений и модели наблюдений, а Нк+1 — матрица т х и (т < и). Введем функции, называемые соответственно невязкой модели процессов и невязкой модели измерений:

гк+1(щ) = Лк+1щ - щк - т/к+1, пк+1(^) = Нк+1щ - фк+1.

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

3(щ; ^1,^2) = <тгк+1^),гк+1М) + <^2Пк+1М,Пк+1М} , (2)

где Ш1 = Щ* > 0, Ш2 = Ш* > 0 — весовые матрпцы п х п и т х т соответственно. Далее везде (.)* означает транспонирование, а < .,. > — скалярное произведение,

2. Решение задачи последовательного вариационного усвоения данных

Лемма 1. Пусть Ш1 = Щ* > 0 Ш2 = Ш2* > 0; тогда минимум целевого функционала (2) задается решением уравнения Эйлера,:

' [Лк+1] * Ш1Лк+1 + [Нк+1] * Ш2Нк+1] Щ = = [Лк+1 ]* Ш1 [щк + т/к+1] + [Нк+1]* Ш2 [фк+1] .

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

Лемма 2. Пусть Ш1 = I, Нк+1 = I,

Ш = 7 ^(Мк+1), Мк+1 € М", Мк+1 = {0 V 1} , а Лк+1 является, трехдиагональной матрицей с элементами

[лк+1], = =1,...,п;

[лк+1] «+1 = -А = 1,...,п - 1;

[лк+1 ]гг-1 = = 2,...,п;

тогда минимум функционала (2) дается решением системы матричных уравнений:

где икЛ1

В1 -1

7мк+1 в

-

С 0 0

и^-Л1 +

В,

-1

тМ^1 В,

А 0' 0 С

икЛ1 -

ц^1

*>к+т/Г1

. 7^1 МкЛ1 .

А, 0 0 С,+1

ик+1 и,+1

+ т/Г1

7^ М^1

2...П - 1;

С„

0

0 А„_1

Цй +

В«

-1

7М"+1 В«

икЛ1

+ т/Г1

(3)

(4)

(5)

к+1

*к+1

В этом случае сопряженная функция ^+1 является вспомогательной и может служить для исследования качества гипотетической модели. Система матричных уравнений (3)-(5), в свою очередь, может быть решена матричной прогонкой.

3. Метод фильтрации Калмана как частный случай последовательного вариационного усвоения данных

Для решения данного класса задач широко используется метод фильтрации Калмана [2]. Предполагается, что состояние системы <ркЛ1 описывается таким случайным вектором, что невязки модели процессов и модели измерений являются нормально распределенными случайными величинами с нулевыми средними и матрицами ковариации и Дк+1 соответственно:

г(<ркЛ1) - N(0, дк+1), ~ N(0, ЯкЛ1).

В этом случае ставится задача минимизации функционала

Ла1шаи(/+*) = Е |/+1 - |

где Е [.] обозначает операцию математического ожидания. Минимум этому функционалу доставляет оценка

/+1 = + КкЛ1 [0кЛ1 - Н, (6)

где

, Г , 1 -1

; (7)

Кк+1 = рк+1 к+1] * Нк+1рк+1 |~Нк+1] * + Дк+1

РкЛ1 = ([Лк+1]-1) |> + ([Лк+1]-1)* ;

Рк = есу(/ - <рк). Матрицу ковариации ошибки прогноза РкЛ1 можно вычислить как

РкЛ1 = есу (/Л1 - (ркЛ1) = = (I - КкЛ1НкЛ1)РкЛ1(/ - КкЛ1НкЛ1)* + КкЛ1дкЛ1 [КкЛ1 ]>

(8) (9)

Лемма 3. Оценка, получаемая по методу Калмана (6)-(10), и оценка, получаемая как минимум функционала вида (2)

3

щ; (Лк+1Рк+1 [Лк+1]*)-1 , [Як+1]

-1

совпадают.

В алгоритмах типа фильтрации Калмана требуется вычисление матриц ковариации ошибки прогноза Рк+1. Хотя это и позволяет оценить качество прогноза, но для ее нахождения требуется обращать полные матрицы, что в случае моделей высокой размерности не всегда возможно, особенно в реальном времени.

4. Расщепленная задача усвоения данных

Пусть модель системы, описываемая уравнением (1), может быть приближена расщепленной системой [3]:

Ак+1фкЦ =(рк + Т^+1.

Кк+1фк+1 =фк+ь.

В этом случае определим решение расщепленной задачи усвоения как последовательное решение задач минимизации следующих функционалов:

1, Ш2) = (т [лк+1(/^ - <рк - , [лк+1(/^ - <рк - ) + (11)

3(щк+1; Ш1,Ш2) = (Ш

Лк+1 Щк+1

Ф

Лк+1 Щк+1

к+±

+

+ <^2П к+1( Щ к+1),п к+1( Щк+1)}.

(12)

Аналогично можно определить расщепленную задачу усвоения данных и для других схем расщепления.

5. Тестовый пример

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

д д д

— (х,у,г) + и(х,у) — (х,у,г) + у(х,у) — (х,у,г) - = /(ж, у, г); (13)

(ж, у, г) € О = [0, Рж] х [0, ¿у] х [0, Т]; (14)

щ(дО, г) = 0, щ(О, 0) = 0, (15)

где /(ж, у, г) задается как точечный источник, расположенный в точке (¿ж/4, ¿у/4) мощности /р и действующий на промежутке времени [0, Т/60]. Задача (13)—(15) аппроксимировалась на равномерной сетке с помощью схемы направленных разностей. Расщепление проводилось покоординатно.

Далее предполагалось, что /(ж, у, £) неизвестна и оценку требуется построить на основе данных измерения функции состояния системы ^ = <р (ж(г), у(г), г) в точке, координаты которой изменяются во времени по закону

ж (г) = Ьж/4, у (г) = £у - г, г < ¿1, у (г) = (г - ¿1), ¿1 < г < г2,

где — скорость движения измерительной системы, В качестве гипотетической модели процесса уз в алгоритме усвоения данных использовался дискретный расщепленный аналог краевой задачи (13)—(15) с /(ж, у, г) = 0,

Параметры модели:

Ьх = 1000 км, Ьу = 1000 км

Т = 20 ч

х у

число точек по времени 300 у = 100

модуль вектора (и, у) = 10 м/с 7 = 100 /р = 3

Бу = 300 км/ч

¿1 = 3.33 ч, ¿2 = 6.66 ч

Рис. 1. Параметры численного эксперимента.

точнее аешение 100 мин оценка решения 100 мин

ЭД®

1001]

2 500

0

1000 | 500 0

1000

3 500

0

1000 2 500 0

0

500 1000 км ■500 мин

0 500 1000 км 000 мин

0

500 1000 км 1200 мин

0

500 1000 км

3 500 0

1000 3 500 0

1000 3 500 0

1000 3 500 0

0 500 Ю00 км ЯШ мин

0 500 1000 км ООО мин

0 500 Ю00 км 1200 мин

0 501] Ю00 км

Рис. 2. Заданное точное решение (слева) и ого оценка последовательной системой усвоения данных (справа).

Для рассмотренной системы ставилась расщепленная задача усвоения. Веса в соответствующих функционалах (11), (12) выбирались как W1 = /, W2 = 7/. Числовые значения параметров эксперимента представлены на рис, 1, На рис, 2 представлены результаты расчетов, Отметим, что алгоритму усвоения удалось качественно воспроизвести по одному измерению на каждом шаге по времени, характер поведения наблюдаемой системы с помощью экономичной схемы вычислений с числом операций O(nx * ny) на один шаг. При этом, однако, абсолютные величины точного решения и его оценки различаются. Успех алгоритма усвоения данных в этом примере объясняется тем, что наблюдательный комплекс в начале временного интервала пересек след от источника.

Заключение

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

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

[1J Пененко В.В. Вариационное усвоение данных в реальном времени // Вычисл. технологии. 2005. Т. 10. Спецвыпуск. Ч. 1 С. 9-20.

[2] Kalman R.K. A new approach to linear filtering and prediction problems // Trans. AME, J. Basic Eng. 1960. Vol. 82. P. 34-35.

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

[4] Penenko A.V. Sequential variational data assimilation // Proc. of SPIE 2005. Vol. 6160. P. 61602D 1-9.

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