Динамико-стохастический подход в задаче усвоения данных наблюдений Текст научной статьи по специальности «Математика»

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

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

ДИНАМИКО-СТОХАСТИЧЕСКИЙ ПОДХОД В ЗАДАЧЕ УСВОЕНИЯ ДАННЫХ НАБЛЮДЕНИЙ

Е. Г. Климова

Институт вычислительных технологий СО РАН, Новосибирск, Россия

e-mail: [email protected]

In this report methods and the data assimilation algorithms based on the Kaiman filter theory are proposed. For calculation of the forecast error covariances the proposed methods use the suboptimal algorithms in which the statistical averaging is replaced by time averaging. A technique of the estimation of the model parameters as well as a "bias" of the forecast model in the data assimilation procedure are considered. The proposed numerical algorithms are compared with the Kaiman filter using a simple test model.

Введение

Задача усвоения данных наблюдений является в настоящее время очень популярной при моделировании различных процессов в атмосфере, океане, водоемах и т. д. Проблема совместного учета данных наблюдений и прогностической модели решается на основе общего оптимизационного подхода с привлечением вариационной постановки либо теории оценивания. В первом случае разрабатывается метод усвоения, называемый в литературе 4DVAR, во втором — алгоритм усвоения, основанный на фильтре Калмана. Если модель рассматриваемых процессов линейна, а случайные величины являются гауссовскими, оба эти подхода математически эквивалентны [1, 2]. Практическая реализация как 4DVAR, так и фильтра Калмана в полной постановке для современных прогностических моделей невозможна из-за высокой размерности возникающих при этом систем уравнений и нелинейности прогнозируемых процессов. Поэтому в обоих случаях рассматриваются различные упрощения алгоритмов. Так, в 4DVAR для вычисления градиента функционала используется линеаризованная модель с упрощенной физикой. Кроме того, в этом подходе изначально считалось, что ковариации ошибок прогноза не меняются со временем. Чтобы учесть эту изменчивость, в настоящее время используются ансамблевый подход и сингулярные векторы линеаризованной модели. Основным направлением в реализации фильтра Калмана является ансамблевый подход [3, 4]. Следует отметить, что этот подход не что иное, как применение метода Монте-Карло для вычисления ковариационных матриц ошибок оценивания. Вторым подходом к практической реализации алгоритма фильтра Калмана является использование субоптимальных алгоритмов, в которых для вычисления ковариаций ошибок прогноза используются упрощенные модели [5 - 8].

Если предположить, что случайные ошибки прогноза обладают свойством эргодичности [9], можно рассмотреть алгоритм, альтернативный ансамблевому фильтру Калмана,

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

при котором вероятностное осреднение заменяется на осреднение по времени. Такой подход предложен в работе [10]. Кроме того, он содержится в работах A.A. Красовского [11]. Алгоритм легко реализуется, однако требуют своего изучения вопросы его применимости в задачах усвоения данных, сходимости алгоритма и связи его с алгоритмом фильтра Калмана. В данной работе рассматривается часть этих вопросов.

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

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

1. Фильтр Калмана

Используемый в работе алгоритм усвоения данных основан на теории оптимальной фильтрации Калмана. Постановка задачи оптимальной фильтрации подробно изложена в [1], описание алгоритма фильтра Калмана можно найти в [2, 5-7].

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

xk = Ak-iXk-i,

где xk — n-вектор прогнозируемых значений в момент времени tk.

Пусть "истинное" значение в момент времени tk описывается векторным уравнением

xk = Ak-ixk-i + £fe_ i, (1) где — случайный вектор "шумов" модели, удовлетворяющий условиям

Eel = 0,

E (4)(£Í)T = Qk öki. Вектор данных наблюдений представим в виде

yo = Mk xk+eo,

где yk — m-вектор наблюдений в момент времени tk; Mk — (nx т)-матрица; £k — случайный m-вектор ошибок наблюдений, удовлетворяющий условиям

E£k = 0,

E (£k)(£k)T = Rk Ski, E (£k)(^)T = 0.

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

хк = Ак-1хк-1, = А-1рк-1АТ-1+Як-1, Кк = Р[ МТ (ЫкР[ МТ + Ек )-1, ра = (1 - КкМк )р[, ха = хк+Кк (у0 - Мкхк),

Рк = Е (хк - хкс )(хк — хкс)Т, Рк = Е (хк — хкс )(хк — х\)Т •

£

Здесь хк — прогноз по модели; х^-1 — вектор проанализированных значений в момент к

времени £ кРк, Рк — матрицы ковариации ошибок прогноза и анализа соответственно; К

оценка анализа совпадает с прогнозом по модели, а матрица ковариаций ошибок анализа равна матрице ковариаций ошибок прогноза.

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

2. Субоптимальный алгоритм, основанный на фильтре Калмана

В предлагаемом в данной работе алгоритме наряду с уравнением для искомой величины рассматривается также и уравнение для ее ошибки:

Дх{ = А к-1Дх{_1 + 4 >

Дх£ = Дх{ - Кк(у0 - Мкхк)•

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

1

N

Ркк = Ахк(Ахк)т = ^Лхг(Ах,)т. (2)

г=1

То есть в данном алгоритме при вычислении ковариационных матриц также использовалось предположение о том, что вероятностное осреднение можно заменить осреднением по времени [9, 10].

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

В работе рассматривался вариант субоптимального алгоритма в применении к задаче оценки систематической ошибки модели. Зададим "истинное" поле в виде

х! = + Ск-гЯк-1 + 4-х,

ок = 0к-1,

где Ак-1 — оператор модели; хк — "истинное" значение в момент времени £к; ок-1 _ "истинная" систематическая ошибка модели в момент времени ¿к-1; матрица Ок-1 в общем случае неквадратная; ошибки е|-1 обладают свойствами

Ее! = 0,

Е (ек)(е|)т = ^ 5к1.

к

хк

хк = Ахк-1 + Ск-1°к-1,

Ахк = АДх+-1Ск-1Аок-1,

к _ к °к = °к-1,

А Ок = А Ок-1-

£

Здесь ок-1 _ прогноз систематической ошибки модели в момент времени £к-1. Параллельно с полями хк, ок вычисляются значения Ахк и Аок, оде Ахк = хк — х{ — отклонение оценки прогнозируемого значения от "истинного", а Аок = Ок — Ок — отклонение оценки систематической ошибки модели от "истинной" систематической ошибки.

Алгоритм оценки параметра ок-1 состоит в том, что вместо исходно го вектора х в

х

этом случае шаг анализа алгоритма фильтра Калмана запишется в виде

хк = хк + Р4,к МТ (Мк Р/ МТ + Як )-1 (уО — Мк хк),

оа = ок + к МкТ (Мк Рк Мкт + Як )-1 (уО — Мк хк),

Ахк = Ахк — р4,кМГ (Мк Р/ МТ + Як )-1(УО — Мк хк),

А ок = А ок — <,к МТ (Мк Ркк МТ + Як )-1(УО — Мк хк),

кк где РХ к — матрица ковариации ошибок прогноза; РХ? к — матрица кроссковариаций ошибок прогноза и оценки систематической ошибки модели, вычисляемые с помощью следующих соотношений:

__1 N

Ри = Лхк(АхкГ = —— ]Г.Ахг(Ахг)т, (3)

— 1 ¿=1

__1 N

Р^М = = Е (4)

В формулах (3), (4) черта сверху означает осреднение по времени.

Представим уравнения для прогнозируемой величины и систематической ошибки в виде (индексы для удобства опущены)

xn+1 = Axn + qn,

qn+1 = qn

Рассмотрим расширенный вектор состояния

и матрицу

A =

q

A 0 0 I

Тогда система перепишется следующим образом:

_ .

Применение алгоритма фильтра Калмана к этой системе на шаге прогноза ковариаций имеет вид

рп+1 _ Дрп ДТ I рп дТ I Дрп I рп

хх хх ' дх ' х^ '

рп+1 _ рп дт I рп

рп+1 _ р п

дд _ дд.

Изменение ковариационных матриц после процедуры анализа алгоритма фильтра Калмана описывается следующими формулами:

р1. _ Р1 - Р1 МТЯ-1 МР1 ,

Pqa = Pf - Pf MTR-1MPf ,

4qx qx ^ qx v ixi

Pa = Pf - Pf MT R-1MPf .

qq qq qx xq

3. Численные эксперименты

Численные эксперименты проводились для простого одномерного уравнения адвекции

dx U дх dt a sin 90 ОХ ' x(0, Х) = F0 sin Х,

где U — константа. Уравнение рассматривалось для круга широты 9 = 45°, брались периодические граничные условия, шаг сетки ДХ = 2.5°, шаг по времени At = 1 ч. Решение искалось с помощью лагранжева подхода. В качестве начального значения бралось Хо(Х) = Fo sin(X), Fo = 10.

Сравнение предлагаемых алгоритмов с фильтром Калмана проводилось с помощью

x

лось, что каждые 6 ч есть данные в подобласти, представляющей собой полосу в 36 узлов сетки. Как это принято в численных экспериментах с модельными данными, "истинное" x

z

значений xl0 бралось x0(Л) = х0(Л) + £0, оде Со = OfN(0,1), Of = 2, N(0,1) — случайная величина, распределенная по нормальному закону с математическим ожиданием, равным нулю, и дисперсией, равной единице. Данные наблюдений задавались наложением на "истинное" значение случайной ошибки наблюдений £0 = o0N(0,1), о0 = 1.

В численных экспериментах проводилось усвоение модельных данных с помощью субоптимальных алгоритмов. В первом эксперименте делалась оценка поля с помощью алгоритма фильтра Калмана и двух субоптимальных алгоритмов. В первом алгоритме ковари-ации ошибок прогноза вычислялись по формуле (2), во втором алгоритме в этой формуле осреднение по времени не проводилось. Результаты эксперимента приведены на рис. 1, где FK, SOI и S02 — это значения среднеквадратической ошибки при применении фильтра Калмана, первого и второго алгоритмов соответственно. Как видно из рисунка, предлагаемые субоптимальные алгоритмы дают поведение ошибки прогноза, близкое к алгоритму фильтра Калмана.

Второй эксперимент проводился для оценки в процедуре усвоения параметра q. В начальный момент искомый параметр q полагался равным нулю, "истинное" значение q задавалось равным 0.01. На рис. 2 приводится значение среднеквадратической ошибки про-

q

2

1.8

1.6

1.4

1.2

0.8

0.6

0.4

0.2

Рис. 1. Среднеквадратическая ошибка прогноза с усвоением.

2

1,3

1,6 1,4 1,2

0,3

0,6 -

мм 11 I I 11 I I 11 11 I I I

т- "в1 Г— О П (О СО п п ^ ^

Время, ч

Рис. 2. Среднеквадратическая ошибка прогноза с усвоением и оценкой параметра, q(0) = 0.

0,0108 0,0106 0,0104 0,0102 0,01

0,0098 -I \

0,0096 -

0,0094 - .......................

0,0092 -0,009 -

0,0088 In 111 11 I 11 11 I 11 111 111 11 i 11 111 111 11 i 11 11 i 111 11 in

T-moionf— T-inoirof— T-in

т-т-счсчсчтп^-^-

Время, ч

Рис. 3. Среднеквадратическая ошибка оценки параметра q.

На рис. 3 представлено поведение среднеквадратпческой ошибки параметра q. Как видно из рисунков, субоптимальный алгоритм и фильтр Калмана дают близкое поведение

q

q

высока, что объясняется отсутствием какой-либо априорной информации о нем в начальный момент времени.

Рассмотренная здесь методика применялась в задаче переноса и диффузии пассивной примеси для оценки эмиссии примеси и систематической ошибки модели. Численные эксперименты проводились совместно с Н.В. Килановой, результаты, полученные в их ходе, опубликованы в работах [14-16].

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

[1] Jazwinski А.Н. Stochastic processes and filtering theory. N. Y.: Acad. Press, 1970. 377 p.

[2] Ghil M., Malanotte-Rizzolli P. Data assimilation in meteorology and oceanography // Advances in Geophysics. 1991. Vol. 33.

[3] Houtekamer P.L., Mitchell H.L. Data assimilation using an ensemble Kalman filter technique // Monthly Weather Review. 1998. Vol. 126. P. 796-811.

[4] Houtekamer P.L., Mitchell H.L. Ensemble Kalman filtering // Quarterly J. of the Royal Meteorological Society. 2005. Vol. 131. P. 1-23.

[5] Климова Е.Г. Упрощенные модели для расчета ковариационных матриц в алгоритме фильтра Калмана // Метеорология и гидрология. 2000. № 6. С. 18-30.

[6] Климова Е.Г. Модель для расчета ковариаций ошибок прогноза в алгоритме фильтра Калмана, основанная на полных уравнениях // Метеорология и гидрология. 2001. № 11. С. 11-21.

[7] Климова Е.Г. Численные эксперименты по усвоению метеорологических данных с помощью субоптимального фильтра Калмана // Метеорология и гидрология. 2003. № 10. С. 54-67.

[8] Todling R., Cohn S. Suboptimal shemes for atmospheric data assimilation based on the Kalman filter // Mon. ma. Rev. 1996. Vol. 124. P. 2530-2557.

[9] Яглом A.M. Корреляционная теория стационарных случайных функций. Л.: Гидрометеоиз-дат, 1981. 280 с.

[10] Сонечкин Д.М. Динамико-стохастический подход к задаче объективного анализа данных разнородных метеорологических наблюдений // Тр. Гидрометцентра СССР. 1976. Вып. 181. С. 54-76.

[11] Справочник по теории автоматического управления / Под ред. A.A. Красовского. М.: Наука, 1987. 711 с.

[12] Dee D.P., Da Silva A.M. Data assimilation in the presence of forecast bias // Quarterly J. of the Royal Meteorological Society. 1998. Vol. 124. P. 269-295.

[13] Балакришнан A.B. Теория фильтрации Калмана. M.: Мир, 1988. 168 с.

[14] Климова Е.Г., Киланова Н.В. Усвоение данных наблюдения в задаче совместного оценивания концентрации и эмиссии пассивной примеси // Матер. V Междунар. сими. "Контроль и реабилитация окружающей среды". Томск, 2006. С. 113-115.

[15] Климова Е.Г., Киланова Н.В. Численные эксперименты по оценке эмиссии метана на основе системы усвоения данных о пассивной примеси в атмосфере Северного полушария / / Оптика атмосферы и океана. 2006. № 11. (В печати).

[16] Киланова Н.В., Климова Е.Г. Численные эксперименты по оценке систематической ошибки модели в задаче усвоения данных о концентрации пассивной примеси / / Вычисл. технологии. 2006. Т. 11, № 5. С. 32-40.

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