Численные эксперименты по оценке систематической ошибки модели в задаче усвоения данных о концентрации пассивной примеси Текст научной статьи по специальности «Математика»

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

Том 11, № 5, 2006

ЧИСЛЕННЫЕ ЭКСПЕРИМЕНТЫ ПО ОЦЕНКЕ СИСТЕМАТИЧЕСКОЙ ОШИБКИ МОДЕЛИ В

ЗАДАЧЕ УСВОЕНИЯ ДАННЫХ О КОНЦЕНТРАЦИИ ПАССИВНОЙ ПРИМЕСИ*

Н. В. КИЛАНОВА, Е. Г. КЛИМОВА Институт вычислительных технологий СО РАН, Новосибирск, Россия e-mail: kilanova@mail.ru, klimova@ict.nsc.ru

An error estimation algorithm for a problem of advection and diffusion of a passive pollution in data assimilation procedure is proposed. Results of the experiments have shown the applicability of the proposed algorithm to the problem of a passive pollution data assimilation. Properties of the algorithm are tested with the help of numerical experiments with the model data.

Введение

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

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

Прогноз полей концентрации дается с помощью полулагранжевой модели переноса и диффузии пассивной примеси для Северного полушария. Для реализации алгоритма

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 04-05-64481).

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

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

1. Модель переноса и диффузии пассивной примеси для Северного полушария

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

ds udsvds ds

+---^т + + (w - Wg) —

dt a cos р dX ao<^ dz

d (1 ds\ 1 d / ds\ d / ds\ „

+ -"^T" cos р— + — kz— + f,

a2 cos2 рдЛ \ dX) a2 cos р др \ dp J dz \ dz J

s(X,p,z,t) = s(X + 2n,p,z,t),

s = sG при p = 0,

ds

v— = a (s — s0) при z = b + h,

dz

ds

v— = 0 при z = H.

dz

Здесь X — долгота; p — широта; z — высота над уровнем моря; s — концентрация примеси; V = (u,v,w — wg) — вектор скорости ветра с компонентами в направлениях X, p,z соответственно; wg — скорость гравитационного оседания; a — средний радиус Земли; k\,k2 — коэффициенты турбулентного обмена в горизонтальном и kz — вертикальном направлениях; f — функция, описывающая размещение и мощность источников; b(X, р) — функция, описывающая рельеф подстилающей поверхности; h — высота приземного слоя и H — верхняя граница расчетной области; параметр a характеризует взаимодействие атмосферы с подстилающей поверхностью [5].

Задача рассматривается в области

Dt = G х [0,T],

где G = S х [b + h, H]; S = {0 < X < 2п; 0 < р < п/2}.

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

Подробнее описание модели и численных экспериментов с системой усвоения для этой модели приведено в [2 - 4].

1

2. Усвоение данных наблюдений в задаче переноса и диффузии пассивной примеси

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

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

хк = Лк-хХк-г. (1)

Здесь Хк — п-вектор прогнозируемых значений концентрации примеси в момент времени ¿к, его размерность равна количеству узлов сетки: п = п\п^пх, где п\ — количество узлов сетки по долготе, п^ — по широте, п — по высоте.

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

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

Ее1 = 0, Е (£* )(£[)т = Qk¿k^

Вектор данных наблюдений представим в виде

ук = Мк хк + а (3)

где ук — т-вектор наблюдений в момент времени ¿к; Мк — (пхт)-матрица; £к — случайный т-вектор ошибок наблюдений, удовлетворяющий условиям

Е& = 0, Е (£к)(£к)Т = ^ ¿н, Е (4к)(£[)Т = 0.

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

хк = Лк_1х2_1; (4)

Р/ = Лк_1 Рк_1ЛТ_1 + Qk_l; (5)

Кк = Р/ МТ (Мк Р/ МТ + Як )-1; (6)

Рк = (/ - Кк Мк )Р/; (7)

ха = хк + Кк (у0 - Мк хк); (8)

Р/ = Е (хк — хк)(х к — хк)Т, Ра = Е (хк — хк )(хк — хк)Т •

/

Здесь х к — прогноз концентрации по модели; хк_1 — вектор проанализированных значений в момент времени £ к _1, т. е. искомая нами оценка, которая получена на (к — 1)-м шаге по

времени; оценку x^-1 называют оценкой анализа; Pf, Pja — матрицы ковариации ошибок прогноза и анализа соответственно; K — весовая матрица.

Таким образом, алгоритм фильтра Калмана состоит из шага прогноза, когда по прогностической модели вычисляется предварительная оценка концентрации по формуле (4) и рассчитываются матрицы ковариации ошибок прогноза по (5). Затем вычисляется весовая матрица K по формуле (6). Далее идет шаг анализа, который заключается в получении непосредственно искомой оценки анализа по (8) и матрицы ковариации ошибок анализа по формуле (7).

В моменты времени, когда наблюдения недоступны, считаем, что оценка анализа совпадает с прогнозом по модели, а матрица ковариации ошибок анализа равна матрице ковариации ошибок прогноза.

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

3. Численные эксперименты с субоптимальным алгоритмом фильтра Калмана

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

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

хк = м- + Gfc-iqt -1+£\-1,

Ок = Qk-1,

где A — оператор модели; х\ — "истинная" концентрация метана в момент времени tk; qk-i — "истинная" систематическая ошибка модели в момент времени tk-1; матрица Gk-1 в общем случае неквадратная; ошибки е\-1 обладают свойствами

Eel = 0; E (е\ )(e[)T = Qköki.

Прогноз концентрации по модели xj, (предварительная оценка поля концентрации) задавался следующим образом:

xk = Axk-1 + Gk-1 qf-1; (9)

Axk = AAxk-1 + Gk-1ÄQk-1; (10)

Qfk = QU (11)

Aqk = AQk-1. (12)

Здесь — прогноз систематической ошибки модели в момент времени ¿к_1. Параллельно с полями хк, вычислялись значения Ахк и А^к, где Ахк = х| — хк — отклонение оценки концентрации от "истинного" значения концентрации, а А^к = вк — вк — отклонение оценки систематической ошибки модели от "истинной" систематической ошибки.

В такой постановке задачи правую часть Ск_1^к_1 можно рассматривать как эмиссию. Результаты численных экспериментов по оценке эмиссии в задаче переноса и диффузии пассивной примеси на примере атмосферного газа метана приведены в [3, 4].

Алгоритм оценки параметра ^к_1 состоит в том, что вместо исходного вектора х в алгоритме фильтра Калмана оценивается расширенный вектор состояния (х,в) [16]. В этом случае шаг анализа алгоритма фильтра Калмана запишется в виде

хк = хк + Р^к МкТ (Мк Р/ МкТ + Як )_1(ук — Мк хк); (13)

вк = в/ + Р4к МТ (Мк Р/ МТ + Як )_1(ук — Мк хк); (14)

Ахк = Ахк — Р^к МТ (МкР/ МТ + Як )_1(ук — Мкхк); (15)

Авк = Авк — / МТ (МкР/ МТ + Як )_1(ук — Мкхк), (16)

где Р/х к — матрица ковариации ошибок прогноза оценки концентрации; Р/ к — матрица кроссковариаций ошибок прогноза оценки концентрации и оценки систематической ошибки модели.

В данном алгоритме при вычислении ковариационных матриц использовалось предположение о том, что вероятностное осреднение можно заменить осреднением по времени [17, 18]. Таким образом, в формулах (17), (18) черта сверху обозначает осреднение по времени:

_ 1 *

Р/х, к = Ах к (Ахк )Т = у Е Ахг(Ахг)Т; (17)

г=1

N

/ = Ахк(Авк)Т = ^ Е Ахг(А9г)Т• (18)

г=1

В численных экспериментах расчеты проводились на сетке 2.5° х2.5° по горизонтали и пятнадцати уровнях по вертикали, шаг по времени А£ = 1 5 мин. Были использованы реальные данные о полях скорости ветра, температуры и давления и модельные данные о концентрации метана [19]. Расчеты проводились для временного периода 48 ч с усвоением каждые 12 ч.

В первом эксперименте в моменты наблюдений моделировались данные наблюдений по формуле (3) таким образом, что в каждый момент усвоения задавалась информация в полосе. На рис. 1 приведено распределение данных наблюдений. При этом каждая вертикальная полоса соответствует области, в которой задаются данные в 12, 24, 36 и 48 часов.

Положим, что в начальный момент времени истинное значение систематической ошибки постоянно для всей области решения задачи. В эксперименте задавалось вО = 0.0009. Кроме того, задано отклонение оценки параметра от истинного значения Адо = 0.0009 х N(0,1 ) х 0.1, где N(0,1 ) — нормально распределенная случайная величина с нулевым математическим ожиданием и дисперсией, равной 1.

Рис. 1. Моделируемое распределение данных наблюдений.

Рис. 2. Среднеквадратическая ошибка: параметра ^, ppm (а) и оценки концентрации, ppm (б). Эксперимент 1.

Результаты эксперимента приведены на рис. 2, а и б. Из рис. 2, а видно, что заданная в начальный момент времени величина Ад0 = д0 — д^ убывает в процессе усвоения, на рис. 2, б отражено убывание среднеквадратической ошибки концентрации.

Во втором численном эксперименте данные наблюдений моделировались в точках реальных спутниковых наблюдений. Координаты точек наблюдений и значения концентраций в них размещены на сайте [20]. На рис. 3 представлено распределение данных наблюдений, которые участвовали в усвоении в момент времени I = 12 ч. По вертикали данные моделировались на четырех верхних уровнях, а именно на высотах 13 660, 16 250, 21 374 и 25 290 м.

Данные наблюдений моделировались таким образом, что невязка в точках наблюдений составляла 0.001. По формулам (9)-(12) определяли предварительную оценку концентрации и систематической ошибки. В моменты усвоения по (13)-(16) получена оценка анализа концентрации и систематической ошибки модели. Результаты этого эксперимента приведены на рис. 4.

Рис. 3. Распределение данных наблюдений в момент времени £ =12 ч.

в

Рис. 4. Среднеквадратическая ошибка: параметра дк, ррт (а); оценки концентрации, ррт (б); оценка значения систематической ошибки модели дк, ррт (в). Эксперимент 2.

Как видно из рис. 4, а и б, среднеквадратические ошибки оцениваемого параметра qk и среднеквадратическая ошибка концентрации убывают после каждого усвоения. Следует отметить, что убывание среднеквадратической ошибки концентрации в эксперименте 1 (см. рис. 2, б) регулярнее и больше по значению, чем в эксперименте 2 (см. рис. 4, б). Это естественным образом объясняется тем фактом, что в каждом усвоении в первом эксперименте участвовало значительно большее количество точек наблюдений, чем в эксперименте 2. Из рис. 4, в видно, что восстановление значения систематической ошибки производится с каждым усвоением, при этом в начальный момент времени значение систематической ошибки равно нулю. В эксперименте с модельными данными показано, что при усвоении происходит убывание ошибки как самой концентрации, так и оцениваемого параметра.

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

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

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

[2] Климова Е.Г., КиллновА Н.В. Усвоение данных наблюдений в задаче переноса и диффузии пассивной примеси // География и природные ресурсы. 2004 (Спецвыпуск): Тр. Меж-дунар. конф. ENVIROMIS-2004. Новосибирск, 2004. C. 175-180.

[3] КилАновА Н.В., Климова Е.Г. Оценка полей концентрации метана над Северным полушарием по данным измерений и модели переноса и диффузии пассивной примеси // Вы-числ. технологии. 2005 (Спецвыпуск): Тр. Междунар. конф. CITES-2005. Новосибирск, 2005. C. 132-137.

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

[5] Марчук Г.И., Алоян А.Е. Глобальный перенос примеси в атмосфере // Изв. РАН. Физика атмосферы и океана. 1995. Т. 31, № 5. C. 597-606.

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

[7] Толстых М.А. Полулагранжева модель атмосферы с высоким пространственным разрешением для численного прогноза погоды // Метеорология и гидрология. 2001. № 4. C. 5-16.

[8] Bermejo R., Staniforth A. The conversion of semi-Lagrangian advection scheme to quasimonotone scheme // Mon. Wea. Rev. 1992. Vol. 120. P. 2622-2632.

[9] Jazwinski A.H. Stochastic processes and filtering theory. N.Y.: Acad. Press, 1970. 377 p.

[10] Menard R., Cohn S.E., Chang L.-P., Lyster P.M. Assimilation of stratospheric chemical tracer observations using a Kalman filter. Pt 1: Formulation // Mon. Wea. Rev. 2000. Vol. 128. P. 2654-2671.

[11] Menard R., Chang L.-P. Assimilation of stratospheric chemical tracer observations using a Kalman filter. Pt 2: %2-validated results and analysis of variance and correlation dynamics // Mon. Wea. Rev. 2000. Vol. 128. P. 2672-2686.

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

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

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

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

[16] Балакришнан А.В. Теория фильтрации Калмана. М.: Мир, 1988. 168 с.

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

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

[19] КрупчАтников В.Н., Крылова А.И. Численное моделирование распределения метана по данным наблюдений на поверхности Земли // Оптика атмосферы и океана. 2000. Т. 13, № 6-7. C. 622-626.

[20] http://haloedata.larc.nasa.gov

Поступила в редакцию 1 августа 2006 г., в переработанном виде — 8 августа 2006 г.