Содержание к диссертации
Введение
Глава 1. Полулагранжева модель переноса пассивной примеси для Северного полушария 16
1.1. Введение 16
1.2. Постановка задачи. Уравнение переноса пассивной примеси. Граничные условия 18
1.3. Квазимоиотонная полулагранжева схема решения уравнения переноса пассивной примеси 20
1.3.1. Решение уравнения переноса 20
1.3.2. Решение уравнений диффузии 24
1.4. Численные эксперименты с моделью 27
1.5. Основные результаты главы 1 42
Глава 2. Алгоритм усвоения данных о пассивной примеси, основанный на теории фильтра Калмана 44
2.1. Введение 44
2.2. Дискретный алгоритм фильтра Калмана 45
2.3. Алгоритмы усвоения данных, основанные на теории фильтра Калмана 50
2.4. Система усвоения данных наблюдений, основанная на модели переноса и диффузии пассивной примеси 52
2.4.1. Блок-схема системы усвоения данных наблюдений о концентрации 52
2.4.2.Особенности этапов реализации субоптимальных алгоритмов усвоения данных о концентрации пассивной примеси 55
2.4.3. Программная реализация су б оптимальных алгоритмов усвоения данных о концентрации пассивной примеси .56
2.5, Численные эксперименты по оценке полей концентрации пассивной примеси 60
2.6. Основные результаты главы 2 71
Глава 3. Методика оценки эмиссии пассивной примеси и систематической ошибки модели, основанная на теории фильтра Калмана 73
3.1. Введение 73
3.2. Алгоритм совместной оценки концентрации и эмиссии пассивной примеси, основанный на алгоритме фильтра Калмана 74
3.3. Численные эксперименты по оценке эмиссии пассивной примеси 78
3.4. Алгоритм совместной оценки концентрации и систематической ошибки модели, основанный на фильтре Калмана 88
3.5. Численные эксперименты по оценке систематической ошибки модели 91
3.6. Основные результаты главы 3 98
Заключение 99
Литература
- Постановка задачи. Уравнение переноса пассивной примеси. Граничные условия
- Решение уравнений диффузии
- Алгоритмы усвоения данных, основанные на теории фильтра Калмана
- Алгоритм совместной оценки концентрации и эмиссии пассивной примеси, основанный на алгоритме фильтра Калмана
Введение к работе
Актуальность.
В последней четверти XX века наблюдается резкое потепление климата Земли. Данное явление имеет частично естественный природный характер, отчасти обусловлено все увеличивающимся антропогенным влиянием. Одной из причин изменения климата планеты является увеличение концентрации парниковых газов в атмосфере. Поэтому актуальной на сегодняшний день является задача оценки концентрации этих газов. Ведущая роль в решении этой задачи отведена математическому моделированию. В настоящее время проблему восстановления пространственно-временного распределения полей газовых примесей в атмосфере принято решать на основе алгоритмов совместного учета прогностических моделей и данных наблюдений о концентрации примесей, называемых алгоритмами усвоения данных.
Для увеличения точности вычисления концентрации примесей возникает необходимость оценки систематической ошибки модели и эмиссии примесей. Поэтому исследование вопросов оценки и разработка алгоритмов для решения подобных задач также являются актуальными в настоящее время.
Цель работы состоит в разработке методики усвоения данных наблюдений, основанной на динамико-стохастическом подходе, для решения проблемы моделирования распространения пассивной примеси в атмосфере с использованием данных натурных наблюдений.
Научная новизна работы.
Предложен новый алгоритм усвоения данных наблюдений в задаче переноса и диффузии пассивной примеси, основанный на предположении об эргодичности случайных полей ошибок прогноза. Проведены численные эксперименты, показывающие, что использование этого алгоритма позволяет улучшить точность оценки полей концентрации. Важным свойством предложенного алгоритма является его экономичность.
Предложена новая методика оценки эмиссии пассивной примеси в процедуре усвоения данных для задачи переноса и диффузии примеси. Проведены численные эксперименты с модельными данными и трехмерной полулагранжевой моделью, которые показали, что предложенная методика позволяет оценивать значения эмиссии по данным наблюдения о концентрации пассивной примеси.
Предложена новая методика оценки систематической ошибки в
процедуре усвоения данных, основанной на алгоритме фильтра Калмана.
Свойства методики проверены с помощью численных экспериментов с
модельными данными.
Достоверность изложенных в работе результатов обеспечивается применением математического аппарата теории оценивания, сравнением полученных результатов тестовых и модельных исследований с результатами других авторов.
Практическая ценность работы
Разработанные численные алгоритмы и реализующие их комплексы программ имеют практическую значимость и могут быть использованы для решения задач усвоения данных наблюдения о концентрации пассивной примеси. Рассмотренные алгоритмы позволяют улучшить точность полей концентрации пассивной примеси и обладают важным свойством экономичности. Предложенные в работе алгоритмы усвоения и оценки систематической ошибки модели могут быть также использованы при моделировании процессов в атмосфере и океане.
Работа выполнена при финансовой поддержке РФФИ (грант ведущих научных школ России № 00-015-98543, грант № 04-05-64481).
На защиту выносятся:
Методика усвоения данных наблюдений для модели переноса и диффузии пассивной примеси с использованием динамико-стохастического подхода.
Алгоритмы усвоения данных наблюдений о концентрации, основанные на теории фильтра Калмана. Алгоритм совместной оценки концентрации и эмиссии пассивной примеси, а также алгоритм совместной оценки концентрации примеси и систематической ошибки модели, обобщающие алгоритм усвоения данных о концентрации.
Система усвоения данных о концентрации пассивной примеси в атмосфере, реализованная в виде комплекса прикладных программ, предназначенного для проведения численных экспериментов по оценки концентрации, а также ряда искомых параметров - эмиссии и систематической ошибки модели.
Апробация работы.
Основные результаты диссертации докладывались на X, XI, XII Рабочих группах "Аэрозоли Сибири" (г. Томск, 2003, 2004, 2005), Международных конференциях по измерениям, моделированию и информационным системам для изучения окружающей среды: ENVTROMIS-2004, 2006 (г. Томск),
Международной конференции "Сопряженные задачи механики, информатики и экологии" (г. Горно-Алтайск, 2004), V Всероссийской конференции молодых ученых по математическому моделированию и информационным технологиям (г. Новосибирск, 2004), Международной конференции и школе молодых ученых по вычислительно-информационным технологиям для наук об окружающей среде: "CITES-2005" (г. Новосибирск, 2005), V Международном симпозиуме "Контроль и реабилитация окружающей среды" (г. Томск, 2006).
Публикации.
Общий объем работ по теме диссертации составляет 15 работ (полный перечень публикаций содержится в диссертации). Из них 6 работ опубликованы в центральной печати, 9 изданы в сборниках тезисов. Работы, в которых содержатся основные результаты диссертации, приведены в конце автореферата (6 работ). Объем этих работ составляет 2.3 печатных листа. Вклад диссертанта 1.15 печатных листов. В том числе, 2 статьи опубликованы в журналах, рекомендованных ВАК для представления основных результатов диссертационной работы (общий объем 0.8, вклад диссертанта 0.4). При выполнении работ, опубликованных совместно с научным руководителем, диссертант принимала участие в обсуждении постановок задач, разработке численных алгоритмов, анализе полученных результатов, подготовке и представлении статей и докладов на конференциях. Ею самостоятельно выполнена программная реализация разработанных численных алгоритмов, проведены расчеты тестовых задач и серия численных экспериментов по оценке концентрации пассивной примеси, эмиссии пассивной примеси и систематической ошибки модели на основе алгоритмов усвоения данных наблюдений.
Структура и объем диссертации.
Постановка задачи. Уравнение переноса пассивной примеси. Граничные условия
Рассматривается модель переноса примеси в сферической системе координат (Х,ср, z), где Х- долгота, р- широта, z- высота, отсчитываемая от уровня моря: при z=b+h, = 0 при z=K
Здесь x=x(X,(p,z,t)- концентрация примеси в точке с координатами (X, (p,z) в момент времени /, u={u,v,w-wg}- вектор скорости ветра с компонентами в направлениях координат X, (р, z, соответственно, wg скорость гравитационного оседания пассивной примеси, ки к2- коэффициенты турбулентного обмена в горизонтальном и ку- в вертикальном направлениях, f=f(X, p,t) функция размещения и мощности источников, а-средний радиус Земли, х0- концентрация примеси на подстилающей поверхности, параметр а характеризует взаимодействие атмосферы с подстилающей поверхностью. Задача рассматривается в области Dt=Gx[0,T], где G=Sx[h+b,H]; 8={0 Л 27г;0 р л/2}. Здесь Ь(Л,ф) функция, описывающая рельеф подстилающей поверхности, h- высота приземного слоя и Я- верхняя граница расчетной области, соответственно [40].
На нижней границе расчетной области поставлено граничное условие, которое позволяет учитывать взаимодействие примеси с подстилающей поверхностью в параметризованном виде [40]. При выводе этого граничного условия сделано предположение, что имеет место аналогия между распределениями примесей и потенциальной температуры [47, 48]. Для нахождения значения XQ В каждом узле сетки записывается уравнение баланса примеси на подстилающей поверхности: где К- число всех наземных источников с координатами (Лі, р,), Q i мощность каждого источника, т.е. количество примеси, поступающей в единицу времени на единичную площадь подстилающей поверхности, Р (i=l,2) - коэффициент, характеризующий взаимодействие примеси с і подстилающей поверхностью (суше соответствует Д =0,01 м/с, а водной поверхности/?-=1 м/с). Учитывая условие (1.2.4), получаем х = , (1.2.7) и J3. -w„+a і vrg где хм = х z=b+h Задача (1,2.1)-(1.2.5) описывает два принципиально различных физических процесса, один из которых является процессом переноса примеси с ее сохранением вдоль траектории, а другой - процессом диффузии. Поэтому для решения сформулированной задачи было естественно использовать метод расщепления по физическим процессам
При этом на каждом временном интервале At=(t ,t J в области G решается задача переноса примеси по траекториям: дх + и дх. v дх-, дх, L + i- + (w-w )—L = o dt acos p дЛ a dtp dz (1.2.8) x, = x при t = t и задача турбулентной диффузии: ( дх V 3dzj дх_ dt 2 2 ял a cos р VA 2 a cos 9 а Эр к cos (з Эр + dz + f (2.1.9) 7+1 7 ж = x, при t = t .
Задача (1.2.8) решается при помощи квазимонотонной полулагранжевой схемы [49] со вторым порядком аппроксимации по пространству и первым порядком по времени [50, 51].
Рассмотрим значение концентрации пассивной примеси в каждой точке области решения (Я, p,z) как скалярную величину х(Я, (p,z,t) (частицу). Тогда уравнение переноса этой величины запишется в виде
Алгоритм решения уравнения (1.3.1) состоит в следующем. Считаем, что в момент времени t = г скаляр находился в начальной точке траектории ы ,<р ,г I внутри некоторой ячейки сетки. Эта траектория такова, что ее конечной точкой в момент времени t = t = t + At является выбранный нами узел сетки U ,<р ,z J (зная конечную точку траектории и скорости движения скаляра, можем найти начальную точку траектории). Тогда, зная поле концентрации в момент времени г , интерполируем значения концентрации из узлов сетки, внутри которой оказался скаляр, в начальную точку траектории \Х ,<р ,z у Найденное таким образом значение концентрации равно значению концентрации в lj+1 j+l j+l конечной точке траектории- узле сетки U ,(р ,z
Итак, в каждый момент времени г в каждом узле сетки аппроксимируем дифференциальное уравнение (1.3.1), описывающее траекторию скаляра x{X,(p,z,t). После аппроксимации имеем разностное уравнение
Решение уравнений диффузии
По формуле (2.3.5) рассчитывалась оценка нормы решения W -rufa ../))2. (1-4.5) где N- количество узлов сетки (в данном эксперименте /У=144 37=5328), І, j- индексы узлов сетки по направлениям 1 и соответственно, x(i,j) значение концентрации в узле сетки с индексами /, / . Изменения нормы решения за двое суток составили 2.5% от заданного в начальный момент значения.
Следующая серия экспериментов проводилась с трехмерной моделью. Рассматривалась модель переноса и диффузии пассивной примеси (1.2.1)-(1.2.5). Численное решение поставленной задачи проводилось на сетке АХ=А(р=2.5 по горизонтали и неравномерной по вертикали, с 15 уровнями, отсчитываемыми в метрах от уровня моря: 50, 300, 600, 900, 1200, 1500, 3040, 4370, 7200, 9220, 11830, 13660, 16250, 21374, 25290. Метод решения описан в пунктах 1.3.1 и 1.3.2.
Расчеты проводились на примере атмосферного газа метана. Это пассивный газ, не вступающий в реакции с другими составляющими атмосферы в течение долгого времени, что позволяет избежать учета химических реакций. Причем, для метана с большой долей точности можно принять, что скорость его гравитационного оседания равна нулю.
Данные о распределении метана были предоставлены нам коллегами из ИВМиМГ СО РАН Крупчатниковым В.Н. и Крыловой А.И. Эти данные получены с помощью глобальной климатической модели, дополненной схемой переноса метана. На основе этой модели проведено усвоение имеющейся информации о концентрации СН4, полученной на сети наземных станций NOAA/CMDL за период 1984-1987 гг. [53].
В численных экспериментах с трехмерной моделью были использованы данные объективного анализа Гидрометцентра России о полях скоростей ветра, давления и влажности за 1-3 августа 2002 года.
Первый эксперимент проводился без учета орографии земной поверхности. Начальное распределение концентрации метана на высоте 300 метров представлено на рис.15.
В расчетах были выбраны константы кг=к7=500 м7с 9 9 горизонтальный коэффициент турбулентного обмена, к3=\ м /с -вертикальный коэффициент турбулентного обмена, функция распределения и мощности источников f=0, скорость гравитационного оседания частиц wg=0, параметр, позволяющий в параметризованном виде учесть влияние подстилающей поверхности, ai=(a/v)= 0.001.
Прогноз поля концентрации метана на 48 часов (ррт). Высота 300 м.
Второй численный эксперимент проводился с учетом орографии земной поверхности. Были взяты данные о высотах земной поверхности относительно уровня моря на всем Земном шаре с шагом по широте и долготе 2 , размещенные на сайте [54]. После этого данные интерполировались на сетку с АЛ=А(р=2.50 по горизонтали. По вертикали значение высоты в каждой точке были заменены на высоту ближайшего расчетного вертикального уровня. Полученное таким образом приближенное поле орографии представлено на рис. 18.
Приближенное поле орографии (метры над уровнем моря).
Начальное распределение концентрации для данного эксперимента получено из предыдущего начального распределения (рис.15), но с учетом орографии. Поле концентрации, заданной в начальный момент времени, на высоте 3040 м приведено на рис. 19. В эксперименте использовались те же поля скорости и направления ветра, что и в предыдущем эксперименте. В расчете выбраны те же константы (горизонтальный и вертикальный коэффициенты турбулентного обмена, скорость гравитационного оседания, функция распределения и мощности источников). Но в отличие от эксперимента без учета орографии, здесь выбран отличный от константы коэффициент ai=(a/v) (см. формулу (1.2.4)), параметризующий влияние подстилающей поверхности. Для вычисления а) был использован программный модуль, предоставленный В.А.Шлычковым.
Алгоритмы усвоения данных, основанные на теории фильтра Калмана
При реализации приведенного выше полного алгоритма фильтра Калмана возникает проблема расчета матриц ковариаций ошибок прогноза, связанная с большой размерностью этих матриц. Допустим, что для отыскания численного решения уравнения (2.2.1) в области решения была построена сетка с N узлами. При решении такой задачи, когда в каждом узле сетки отыскивается одна скалярная величина, например задачи (1.2.1)-(1.2.5), размерность п векторов прогноза и анализа равна N. Если же в каждом узле сетки требуется найти векторную величину, состоящую из р скалярных, то тогда n=N p. Размерность матрицы ковариаций ошибок прогноза РІ, вычисляемой в алгоритме фильтра Калмана по формуле (2.2.9), равна (пхп), что при решении рассматриваемой задачи (1.2.1)-(1.2.5) на сетке, описанной в параграфе 1.4, составляет 79920x79920 (№=144 37 15=79920 - произведение количества узлов по направлениям X, p, z, соответственно). В задачах с большим пространственным разрешением размерность этой матрицы может быть в разы или даже на порядки больше. Понятно, что на современных вычислительных машинах рассчитывать такие матрицы невозможно. Это одна из проблем, возникающих при практической реализации алгоритма фильтра Калмана в задачах усвоения. Потому многие исследователи рассматривают в своих работах различные способы нахождения матрицы ковариаций ошибок прогноза, применяя упрощенные операторы для ее вычисления [17, 18-22], либо используя для этого ансамблевый подход [25,26].
Алгоритмы усвоения данных наблюдений, в которых есть какие-либо неточности в задании или вычислении ковариационных матриц (по сравнению с классическим алгоритмом фильтра Калмана), называются субоптимальными алгоритмами усвоения данных [17,18-22].
В данной работе рассматривалось несколько подходов при расчете матрицы ковариаций ошибок прогноза.
1. Первый подход. Для расчета ковариаций ошибок модели используется упрощенный оператор модели. При этом делается предположение, что ошибки прогноза на разных вертикальных уровнях не коррелируют между собой, тогда матрица ковариаций ошибок прогноза становится блочно-диагональной. В этом случае каждый блок можно рассчитывать по двумерной модели переноса, т.е. на каждом уровне по высоте свой блок (шаг диффузии в расчете опущен).
2. Второй подход основывается на предположении об эргодичности случайных полей ошибок прогноза, то есть при вычислении характеристик этих полей теоретико-вероятностное осреднение можно заменить осреднением по времени [56]. В этом случае рассматривается дополнительное уравнение для вычисления ошибки прогноза по модели. При этом матрица ковариаций вместо (2.2.9) вычисляется по формуле //= ( /= ХЛ Лх/, (2.3.1) где Ax = x, - xi - ошибка прогноза в момент времени t, .
Черта сверху означает осреднение по времени. Такой подход к задаче усвоения был предложен Сонечкиным А.Д. [7]. Кроме того, на таком подходе основаны методики оценивания, рассматриваемые в работах Красовского А.А. и др. [57]
Предложенные в работе субоптимальные алгоритмы усвоения данных о пассивной примеси использовались в задаче оценки концентрации метана в атмосфере Северного полушария. Для реализации этих алгоритмов была разработана система усвоения данных, основанная на модели переноса и диффузии примеси.
В численных экспериментах по некоторой априорной информации давался прогноз по модели на 48 часов с усвоением данных каждые 12 часов. Таким образом, четыре раза реализовывались все шаги алгоритма (шаг прогноза, шаг анализа). Блок-схема системы усвоения, по которой проводились численные эксперименты, представлена на рис.22.
Из рис.22 видно, что на предварительном этапе необходимы начальное поле концентрации примеси, а таюке задание статистических характеристик ошибок прогноза и наблюдений. Для субоптимального алгоритма, основанного на первом подходе из параграфа 2.3, это задание f начальной матрицы ковариаций ошибок прогноза Р: , а также матрицы ковариаций ошибок наблюдения R,, необходимой для вычисления весовой матрицы на шаге анализа (см. (2.2.12)), и матрицы ковариаций ошибок модели Q, . Для алгоритма, основанного на втором подходе из f параграфа 2.3, это задание начального поля ошибки прогноза Axi и также задание матриц ковариаций ошибок наблюдения R, и ковариаций ошибок модели Q, .
Далее следует шаг прогноза по модели. В проведенных экспериментах по усвоению данных о пассивной примеси для получения прогноза поля концентрации использовалась трехмерная модель переноса и диффузии пассивной примеси для Северного полушария в сферической системе координат. Модель подробно описана в параграфе 1,2. Вычислялись поля прогноза концентрации на двое суток. Шаг по времени составил 15 мин. Кроме того, при реализации субоптимального алгоритма, основанного на первом подходе из параграфа 2.3, на шаге прогноза вычислялась матрица ковариаций ошибок прогноза по двумерной модели переноса с шагом по времени один час. А при реализации алгоритма, основанного на втором подходе из параграфа 2.3, вычислялся прогноз ошибки оценки концентрации (см. формулу (2.5.3)) по той же модели, что и сама концентрация, то есть по трехмерной модели переноса и диффузии пассивной примеси.
Алгоритм совместной оценки концентрации и эмиссии пассивной примеси, основанный на алгоритме фильтра Калмана
По методике, описанной в параграфе 3.3, была проведена серия численных экспериментов по оценке систематической ошибки модели.
В первом эксперименте задавалось начальное "истинное" значение систематической ошибки q =0.0009. Кроме того, задавалась начальная ошибка оценки параметра Ад -до N(0,1), где #=0.0009 предполагаемый уровень систематической ошибки, Лґ(0,1) нормально распределенная случайная величина с нулевым математическим ожиданием и дисперсией 1, о =0.1 - предполагаемое стандартное отклонение ошибок прогноза. Предварительная оценка систематической ошибки равна нулю: f д -0. Данные наблюдений моделировались в моменты наблюдений (в 12 часов, 24 часа, 36 часов, 48 часов) с распределением по пространству, приведенным на рис.22.
Аналогично экспериментам, описанным в параграфе 2.3, были использованы данные объективного анализа Гидрометцентра России о полях скорости и направлении ветра, давления и влажности за 1-3 августа 2002 года, а также начальное распределение концентрации метана. Для получения оценки концентрации метана в процессе усвоения данных задавалась начальная ошибка оценки концентрации, такая же, как f и в экспериментах, описанных в параграфе 2.5: Ах ха N(0,l)k , где 0 0 е х=0.6 - среднее значение концентрации по расчетной области, и =0.1 стандартное отклонение ошибок прогноза, Я(0,1)-нормально распределенная случайная величина с нулевым средним и дисперсией 1, к =1.2 - эмпирический коэффициент для моделирования начальной е ошибки оценки концентрации. Матрица ковариации ошибок наблюдений R, задавалась по формуле (2.5.2): R? =а21. к 2 Здесь /- единичная матрица, а =0.01 - предполагаемая дисперсия ошибок наблюдения. Матрица ковариации ошибок модели Q, задана нулевой.
В процессе усвоения были получены оценки концентрации и систематической ошибки модели. На рис.45 и рис.46 приведены результаты расчетов. обусловленный добавлением систематической ошибки модели q при вычислении по формуле (3.4.1). При этом можно заметить общее снижение уровня среднеквадратической ошибки оценки концентрации после четырех шагов анализа.
В следующем численном эксперименте производилась оценка параметра систематической ошибки, когда в точках реальных спутниковых наблюдений [57] моделировались данные наблюдений о концентрации (эксперимент 2). По вертикали данные моделировались на четырех верхних уровнях, а именно, на высотах 13660, 16250, 21374, 25290 метров. Пример пространственного распределения точек наблюдения, в которых моделировались данные наблюдений в момент времени ґ=12 часов, приведен на рис.30. Точки, в которых моделировались данные в момент времени /=24 часа показаны на рис.31. В остальном проводимый численный эксперимент аналогичен предыдущему.
При проведении эксперимента было получены оценки среднеквадратических ошибок концентрации и параметра q . На рис.47 к приведена среднеквадратическая ошибка оценки концентрации. Из рисунка видно, что между шагами анализа наблюдается рост среднеквадратической ошибки (как и в предыдущем эксперименте) и в моменты анализа имеет место сокращение ошибки. На рис.48 приведена среднеквадратическая ошибка оценки параметра q . Этот рисунок к отражает убывание в процессе усвоения данных о концентрации величины среднеквадратической ошибки параметра q .
На рис.49 приведены графики средиеквадратической ошибки концентрации, полученные в результате проведения эксперимента 1, в котором данные наблюдений моделировались с распределением по пространству, показанному на рис.22 (сплошная линия) и в результате проведения эксперимента 2, в котором данные наблюдений моделировались в точках реальных спутниковых наблюдений (пунктирная линия). Сравнение двух графиков позволяет сделать вывод, что существенную роль здесь играет количество данных наблюдений. В первом эксперименте на каждом шаге анализа данные имелись в 2880 узлах сетки, а во втором - в 32 точках наблюдений. Причем, проведение шагов анализа в эксперименте 1 позволяет сократить начальный уровень средиеквадратической ошибки оценки концентрации, а во втором - нет (рис.49).
