CC BY
9
УДК 504.064.37
ЧИСЛЕННОЕ ИССЛЕДОВАНИЕ АЛГОРИТМА ВАРИАЦИОННОГО УСВОЕНИЯ ДАННЫХ МОНИТОРИНГА ХИМИЧЕСКОГО СОСТАВА АТМОСФЕРЫ В УСЛОВИЯХ ГОРОДА НОВОСИБИРСКА
Алексей Владимирович Пененко
Институт вычислительной математики и математической геофизики СО РАН, 630090, Россия, г. Новосибирск, ул. Академика Лаврентьева, 6, кандидат физико-математических наук, старший научный сотрудник; Новосибирский национальный исследовательский государственный университет, 630090, Россия, г. Новосибирск, ул. Пирогова, 2, научный сотрудник, тел. (383)330-61-52, e-mail: a.penenko@yandex.ru
Александр Владимирович Гочаков
Сибирский региональный научно-исследовательский гидрометеорологический институт (ФГБУ СибНИГМИ), 630099, Россия, г. Новосибирск, ул. Советская, 30, старший научный сотрудник, тел. (961)874-33-93, e-mail: gochakov@sibnigmi.ru
Жадыра Саматовна Мукатова
Новосибирский национальный исследовательский государственный университет, 630090, Россия, г. Новосибирск, ул. Пирогова, 2, студент; Институт вычислительной математики и математической геофизики СО РАН, 630090, Россия, г. Новосибирск, ул. Академика Лаврентьева, 6, техник, тел. (903)997-81-89, e-mail: zmukatova@ya.ru
Павел Николаевич Антохин
Институт оптики атмосферы им. В. Е. Зуева СО РАН, 634055, Россия, г. Томск, пл. Академика Зуева, 1, научный сотрудник, тел. (3822)49-10-23, e-mail: apn@iao.ru
В работе представлены результаты численного исследования эффективности вариационного алгоритма усвоения данных для модели переноса и трансформации атмосферных примесей. Алгоритм основан на схеме расщепления. На стадии переноса используется прямой алгоритм вариационного усвоения данных, а на стадии трансформации - итерационный. Алгоритм показал свою эффективность в численном эксперименте для города Новосибирска.
Ключевые слова: усвоение данных, перенос и трансформация примесей, вариационный подход, мониторинг загрязнений воздуха, схема расщепления.
NUMERICAL STUDY OF THE ALGORITHM OF VARIATIONAL
DATA ASSIMILATION FOR MONITORING THE CHEMICAL COMPOSITION
OF THE ATMOSPHERE IN THE CITY OF NOVOSIBIRSK
Alexey V. Penenko
Institute of the Computational Mathematics and Mathematical Geophysics SB RAS, 6, Prospect Аkademik Lavrentiev St., Novosibirsk, 630090, Russia, Ph. D., Senior Researcher; Novosibirsk National Research State University, 2, Pirogova St., Novosibirsk, 630090, Russia, Researcher, phone: (383)330-61-52, e-mail: a.penenko@yandex.ru
Alexander V. Gochakov
Siberian Regional Research Hydrometeorological Institute (fgbi, Shinigmi), 30, Sovetskaya St., Novosibirsk, 630099, Russia, Senior Researcher, phone: (961)874-33-93, e-mail: gochakov@sibnigmi.ru
Zhadyra S. Mukatova
Novosibirsk National Research State University, 2, Pirogova St., Novosibirsk, 630090, Russia, Student; Institute of the Computational Mathematics and Mathematical Geophysics SB RAS, 6, Prospect Аkademik Lavrentiev St., Novosibirsk, 630090, Russia, Technician, phone: (903)997-81-89, e-mail: zmukatova@ya.ru
Pavel N. Antokhin
Institute of atmospheric optics. Atmospheric optics SB RAS, 1, Academician Zuev Square, Tomsk, 634055, Russia, Researcher, phone: (3822)49-10-23, e-mail: apn@iao.ru
The results of a numerical study of the efficiency of the variational data assimilation algorithm for the transport and transformation model of atmospheric pollution are presented. The algorithm is based on a splitting scheme. The direct algorithm of variational data assimilation is used at the transport stage, and the iterative algorithm is used at the transformation stage. The algorithm shown its effectiveness in the numerical experiment for the city of Novosibirsk.
Key words: data assimilation, transport and transformation models, variational approach, air pollution monitoring, splitting scheme.
Введение
Методы математического моделирования с оперативным усвоением данных измерений уровня загрязнения, полученных на стационарных или мобильных постах мониторинга позволяют получать информацию о распределении и динамике загрязняющих веществ в атмосфере города. Обзор современных алгоритмов усвоения данных можно найти в [1].
Данная работа обобщает [2] на случай двух пространственных переменных. Современные модели химии атмосферы обычно являются трехмерным по пространству, однако их использование, особенно в контексте обратного моделирования, сопряжено с необходимостью задействования суперкомпьютеров. Для работы на машинах со средними харакетристиками интерес представляют двумерные по пространству модели.
Усвоение одного и того же набора данных измерений в рассматриваемом алгоритме производится квази-независимо на отдельных стадиях схемы расще-пления.Таким образом на каждом шаге модели параллельно решается три группы задач варационного усвоения данных. На каждой стадии схемы расщепления, на ограничениях математической модели, находится условный минимум целевого функционала, совмещающего невязку между измеренными значениями и их смоделированными аналогами, а также некоторый стабилизатор, включающий норму функции неопределенности (управления).
На стадии переноса рассматривается один шаг по времени математической модели. Для такой постановки в [3] рассмотрен и теоретически исследован прямой (без итераций) алгоритм усвоения данных. Численные схемы для реализации усвоения на стадии трансформации были представлены в [4]. При этом схемы для получения прямых и сопряженных задач, необходимых для построения градиентов целевого функционала согласованы в смысле тождества Ла-
гранжа. Целью работы является проверка эффективности алгоритма усвоения данных на правдоподобном сценарии для города Новосибирска.
Методы и материалы
Рассмотрим прямоугольную пространственно-временную область В:=[0, X] х[0, У], (х, г) еОТ := В х [0,Т ], ограниченную дПТ = дОх [0,Т ]. В этой области решается задача о транспорте и трансформации набора химических веществ в атмосфере:
дф дг
+ й1у(иф - №тайф1) + Р1 (ф)ф = П1 (ф) + / 1 + г, (х,г) е0.Т. (1)
^^т = 0, (х,г)еГо»г; ф = (Ф), (хг)еГп; ф(,0) =(Ф,), хеа (2) дп
где ф е ММс - вектор-функция состояния, различные элементы которой ф, I = 1,...,Ыс представляют концентрации рассматриваемых примесей, Ыс - их количество, р(ф)>0 - коэффициенты деструкции рассматриваемых веществ, П (ф)>0 - скорости продукции, и(х, г) е М3 - вектор скорости ветра, Гоиг -части границы области, где скорость ветра направлена наружу из области, Гп -части границы, где скорость ветра направлена внутрь области, /и(х, г) -коэффициент диффузии, /г, (ф0) - априорные значения источников и начальных данных, (фь) - концентрация вещества на границе области, г -
управляющая функция (неопределенность), которая добавляется в жесткую структуру модели для усвоения данных. Операторы продукции и деструкции задаются системой кинетических уравнений Прямая задача состоит в определении ф из (1)-(2) по известным /г,(ф0) ,(фь) , гг. Предположим, что
доступны результаты контактных измерений концентраций {1т }т=1 м
в заданные моменты времени {вт }т=1,...,м в точках {^т }т=1,...,м для заданных
веществ {1т }т=1,...,м . Измерения могут содержать погрешности {^т }т=1м :
4 = ф (Хт ,вт ) + Лт , т = 1,..., М. (3)
т
Ошибка измерений Г!т предполагается ограниченной во (взвешенной)
Эвклидовой норме в пространстве измерений \щм <8^. Мы предполагаем, что все функции и модельные параметры достаточно гладкие, чтобы решения существовали и преобразования имели смысл. Задачу определения ф для г > г * по (1)-(2) и (3) при заданных функциях /г, (ф,) , (фь) и результатах измерений
{1т }т=1 м таких, что 0< вт <г* назовем задачей усвоения данных.
временном интервале
гк -1, гк
Для решения многомерных задач применяется метод расщепления. Мы будем использовать аддитивно усредненные схемы расщепления. Основная идея подхода состоит в том, чтобы усваивать данные локально на отдельных шагах
расщепления. Введем временную сетку щ = |о< г1 <... < гк <... < = . На каждом
приблизим общую модель конвекции-диффузии-реакции (1)-(2) схемой расщепления по физическим процессам и пространственным переменным, порожденной разбиением = 1, / > 0, 2 = {х, у,г}.
Для примера рассмотрим шаг расщепления, соответствующий только одному пространственному измерению х (* = х) для любого I = 1,..., N:
/Т ^/^ + (Х,')е[0'Х]Х (4)
м-^п- =°, (х,%) е Гомг п ({х = 0} и {х = X}) X [гк-1,гк], (5)
дп
тР = т),, (х,г) £ГЙЙ( п ({х = 0} и{X = X})X[гк-1,гк], (6)
тГ{, гк-1)=тХ, гк-1). (7)
В качестве решения задачи усвоения данных рассмотрим минимум целевого функционала
С) ((*"')).(г)) ) = С((«')), - 'к, )2 Мк, + Р±((Г«>)) )2
на ограничениях (4) - (7), аппроксимированных трехдиагональной неявной
разностной схемой, считая
гк -1, гк
одним шагом по времени:
- А ),+с. (*<й)), - ^ =(ф )Г+А'рк+д' (г(й 1>.=1
Здесь с - параметр усвоения, Мк]п1 =1 , если в узле с координатами (хг, у , гк) есть измерение вещества 1, иначе М,п1 =0. Аналогично, 1,п1 равно результату измерения, если в узле с координатами (х., у,, гк) есть измерение вещества 1, иначе оно равно 0. Решение этой задачи оптимизации задается трех-диагональным матричным уравнением [3], которое можно решить прямым методом матричной прогонки.
Условный минимум функционала ф назовем решением задачи усвоения данных на стадии переноса.
В силу жесткости системы дифференциальных уравнений, модель процессов трансформации (» = с)
д »)
у.ф + Р (ф )ф = П1 (ф) + » + , (х, г) е В X [гк-1, гк ],
гк-1) =ф (.,гк-1), х е В,
требуется решать на более тонкой сетке по времени, нежели модель для процессов переноса. Введем на каждом шаге
гк-1, гк
равномерную подсетку
]-0.5 _
< гр <... < г_г = гк | с точками и длиной интервала А г .
[гк-1= 71<...
Для решения системы используется монотонная дискретно-аналитическая схема первого порядка из семейства схем, представленных в [5], аналогичная
части схемы ОББЛ [6]. Для I = 1,..., _с при фр = (ф{»)]и , = г]»:
4Р+1 = ЛР
= фрьр(рр(Трф)) + Ор(Р(Гр,фр))(П(Гр,фр) + ,), р = 1,...,N -1, (8)
р1=ф)!,
(9)
П (Р ) =
0
- РА Г
РА г > г
е'" РАг е[гт1П,гтах], 1 - Р АТ 0< Р АТ < <1
(Р) =
1/Р
РАг > гт
1 - е-РАг -
— Аг РА г е[в1а1п,втШХ ]
РАг АТ
0<РАг <£т. <1
Данные измерений на сетке ю/-05 будем интерпретировать как значение 11 концентраций вещества I в точке гк, при этом функцию управления г]» будем считать постоянной на -1, гк ]. Рассмотрим целевой функционал
N. _ ,2 N
(» )) (ф V) = ((ф<#)) - ] )мкк, + № )2,
(10)
где ф} - решение (8)-(9), м] равно 1, если в узле с координатами (хг, у, гк)
есть измерение вещества 1, иначе оно равно 0, 1кп1 равно результату измерения,
если в узле с координатами (х{, у}, гк) есть измерение вещества 1, иначе оно
равно 0. С помощью сопряженных уравнений можно вычислить градиент (10):
( \\ _г-1 V,)) ((ф»)) ,(г(»)) ) = ^(г]ф+2Р1 (г(»)), 1 = 1,..., _с,
где ^ - решение сопряженной задачи
-1 = ((ф,т>))У, ] = 1,...,N -1,
; Д(^-\/-1)= О(р), р = /(/к-1,фк-1),
Ж(/к-1,фк-1,г) = Лш^ Ь(рг2) + {(,фк-1,г)}
= а/шёь (р] )+{$ (г )}}.
1 =1 N =1
Здесь {$ }/=с1 обозначает матрицу со строками Sг
$ (,ф, г ) = О (р1 )УП 1 (ф,ф)+(ффЬ'(р) + О'(р )(П, (/,ф) + г ))У/> (,ф,фф
Аналогично [7], при минимизации функционала используется регуляризация через проекцию градиентов на вещества, которые соответствуют первичным загрязнителям. После чего задача оптимизации решается методом сопряженных градиентов. Обозначим ее минимум через ф(г).
В качестве оценки решения задачи усвоения данных, примем усреднение на шаге /к
Для подготовки коэффициентов модели переноса примесей в составе алгоритма усвоения данных использовались метеорологические величины, рассчитанные с помощью модели "КБ-СИет [8] в области, ограниченной координатами 54.75° - 55.16° с.ш. и 82.66° - 83.37° в.д. (рис. 1). Параметры сеток в экспериментах с усвоением данных: N = 25 , =720х 5, т = 24 х 3600 . Для таких параметров расчет 24 часов модельного времени на настольном компьютере средних характеристик занимает примерно 2 часа.
• •
Рис. 1. Расположение основных источников (треугольники) и постов мониторинга (круги)
Для описания механизмов трансформации использована модификация схемы из [9], содержащая 22 вещества и 20 химических реакций. Скорости фотохимических реакций зависят от времени суток.
В качестве информации о расположении источников Г использовались данные о городских ТЭЦ и котельных (рис. 1, отмечены треугольниками). Все источники считались одинаковыми с постоянной мощностью по времени и выбросами СО, N0, Б02. С заданным Г решалась прямая задача для генерации
«точного значения» полей концентраций примесей ф . Измерения концентраций С02, £02,03, N02, проводятся с шагом в 2 минуты на 9ти пунктах мониторинга города (рис. 1, отмечены кругами) [10].
Результаты
По отношению к задачам обратного моделирования химические вещества можно разделить на несколько пересекающихся множеств: вещества, которые выбрасываются источниками, вещества концентрации которых измеряются, и «ненаблюдаемые вещества».
На рис. 2 представлен результат усвоения данных, «точное решение» и результат без усвоения в сценарии для города Новосибирска для вещества из группы «ненаблюдаемых» (НСН0) в один из моментов времени.
НСНО 1=11.67
усвоения
[У1: ■
ишг
; • •». •••
Нп) 1 "-г :.:•:• ... .V...
ЩЁ У; Д ' Н^УЛ
■ ; \
(6)
Рис. 2. Сравнение «точного решения» (а) с решением без усвоения (б) и с решением с усвоением данных (в)
На рис. 3 представлены относительные ошибки для вещества из группы «ненаблюдаемых» (НСНО) и вещества, которое одновременно и выбрасывается и «наблюдается» (8О2).
нсно
без усвоения усвоение
302
5 10 б)
без усвоения усвоение
15
20
-I, ч.
Рис. 3. Сравнение относительных ошибок для решений с усвоением данных и без него для различных химических веществ: а) для НСНО, б) для 8О2
Обсуждение
Анализируя рис. 2, можно сделать вывод, что качественно решение с усвоением данных наблюдений лучше соответствует «точному решению», нежели решение без усвоения. Количественные оценки сходимости алгоритмов (рис. 3) подтверждают выводы, полученные качественно при рассмотрении рис. 2. Можно также отметить, что ошибка в «ненаблюдаемом» НСНО начинает возрастать после «восхода солнца» в 5 часов, когда «запускаются» фотохимические реакции. Для «выбрасываемого» и «наблюдаемого» вещества 8О2 связь между «восходом солнца» и поведением ошибки не возникает.
Таким образом, можно заключить, что в рассмотренном численном эксперименте рассмотренный алгоритм усвоения оказался более эффективным, чем алгоритмом без усвоения данных.
Заключение
Таким образом, усвоение данных сети мониторинга позволяет улучшить качество моделирования. Используя схему расщепления, можно построить вычислительно эффективный алгоритм для усвоения данных контактных измерений для моделей конвекции-диффузии-реакции, способный при небольших разрешениях сеток работать в реальном времени на вычислительных машинах средних характеристик, что важно, например, при проведении полевых экспериментов.
Благодарности
Разработка новой модификации прямого алгоритма усвоения данных выполнена в рамках темы госзадания ИВМиМГ СО РАН № 0315-2016-0004, построение сценария для моделирования, адаптация алгоритма к условиям города Новосибирска и численные эксперименты выполнены при поддержке РФФИ и Правительства Новосибирской области в рамках научного проекта № 17-41543309.
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Bocquet M., Elbern H., Eskes H., Hirtl M., Carmichael G.R., Flemming J., Inness A., Pagowski M., Cama, et al. Data assimilation in atmospheric chemistry models: current status and future prospects for coupled chemistry meteorology models // Atmospheric Chemistry and Physics Discussions, -2014. -Vol. 14.- P. 32233-32323.
2. Penenko A., Penenko V., Tsvetova E., Grishina A., Antokhin P. Sequential Variational Data Assimilation Algorithms at the Splitting Stages of a Numerical Atmospheric Chemistry Model // In: Large-Scale Scientific Computing. Springer International Publishing.- 2018.-P. 536-543.
3. Penenko A.V., Penenko V.V., Tsvetova E.A. Sequential data assimilation algorithms for air quality monitoring models based on a weak-constraint variational principle // Numerical Analysis and Applications. -2016. -Vol. 9. -P. 312-325.
4. Penenko A.V. Consistent Numerical Schemes for Solving Nonlinear Inverse Source Problems with Gradient-Type Algorithms and Newton-Kantorovich Methods // Numerical Analysis and Applications. - 2018.-Vol. 11. -P. 73-88.
5. Penenko V.V., Tsvetova E.A., Penenko A.V. Variational approach and Euler's integrating factors for environmental studies // Computers & Mathematics with Applications. - 2014. -Vol. 67. - P. 2240-2256.
6. Hesstvedt E., Hov O., Isaksen I.S.A. Quasi-steady-state approximations in air pollution modeling: Comparison of two numerical schemes for oxidant prediction // International Journal of Chemical Kinetics. - 1978.-Vol. 10.-P. 971-994.
7. Пененко А. В., Гришина А. А. Численное исследование алгоритмов оценки источников для модели трансформации примесей в атмосфере // Интерэкспо ГЕО-Сибирь-2017. XIII Междунар. науч. конгр. : Междунар. науч. конф. «Дистанционные методы зондирования Земли и фотограмметрия, мониторинг окружающей среды, геоэкология» : сб. материалов в 2 т. (Новосибирск, 17-21 апреля 2017 г.). - Новосибирск : СГУГиТ, 2017. Т. 1. - С. 136-140.
8. Skamarock W.C., Klemp J.B., Dudhia J., Gill D.O., Barker D.M., Duda M.G., Huang X.Y., Wang W., Powers J.G. A Description of the Advanced Research WRF Version 3, Boulder, Colorado, USA, 2008.
9. Stockwell W.R. Comment on Simulation of a reacting pollutant puff using an adaptive grid algorithm by R.K. Srivastava et al. // Journal of Geophysical Research. - 2002. - Vol. 107.-P. 4643-4650.
10. Департамент природных ресурсов и охраны окружающей новосибирской области. О состоянии и об охране окружающей среды Новосибирской области в 2015 году, Новосибирск, Государственный доклад 2016. [Online]. http://www.nso.ru/sites/test.new.nso.ru/ wodby_files/files/wiki/2014/01/korrektura_gosdoklad-2015.compressed.pdf.
© А. В. Пененко, А. В. Гочаков, Ж. С. Мукатова, П. Н. Антохин, 2018
