CC BY
26
Вычислительные технологии
Том 11, часть 3, Специальный выпуск, 2006
ДИНАМИКО-ВЕРОЯТНОСТНОЕ МОДЕЛИРОВАНИЕ ВЫБРОСОВ В АТМОСФЕРЕ С ИСПОЛЬЗОВАНИЕМ ПОЛУЛАГРАНЖЕВОЙ МОДЕЛИ ПЕРЕНОСА СУБСТАНЦИИ*
А. В. Протасов Институт вычислительной математики и математической геофизики СО РАН, Новосибирск, Россия e-mail: sasha@ommfao.sscc.ru
For a given local area of size 10x10 of degrees of the Northern hemisphere with use of a method of dynamic probabilistic modeling and six-hour spatial — temporal data of the reanalysis for 1948-2005 a climatic ensemble of spatial — temporal realizations of fields of an meteoelements is constructed. The realizations of the ensemble are solutions of the appropriate numerical model of hydrothermodynamics of the atmosphere and have a correlation matrix close to the correlation matrix estimated by the reanalysis data within the solving of the appropriate variational assimilation problem. Using a definition of a climatic trajectory of atmospheric movements, emission of concrete realizations of trajectories by limits of the climatic trajectory are investigated. These trajectories by are naturally determined by appropriate realizations from the climatic ensemble on the basis of solving of the semi-Lagrangian transport model and the method of detector functions.
Введение
На основе вариационного метода усвоения информации с помощью математической модели и статистических методов моделирования рассматривается метод динамико-вероят-ностного моделирования комплексов многомерных гидрометеорологических полей [1-3]. Реализации моделируемых полей удовлетворяют априорным статистическим свойствам исследуемых полей, заданных в виде ковариационных и взаимно-ковариационных матриц, распределения вероятности и т.д., а также физическим свойствам, присущим решениям соответствующих систем дифференциальных уравнений. По существу, предлагается метод динамико-вероятностного моделирования многомерных гидрометеорологических полей, в общем случае — негауссовских, с оптимальным соответствием статистических и физических свойств.
Гидротермодинамическая модель, заданная системой дифференциальных уравнений, служит интерполятором в пространственно-временные точки рассматриваемой области,
* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 05-05-98000-р^обь^а).
© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2006.
она предназначена также для их взаимного согласования и фильтрации нефизических компонентов, а статистическая модель обеспечивает заданную вероятностную структуру рассматриваемого процесса и определяет ансамбль независимых реализаций полей. Использование вариационного метода усвоения информации позволяет оптимизировать процесс объединения динамических и статистических методов численного моделирования.
Одним из возможных применений предложенного метода является метод численного моделирования климата атмосферы в виде ансамбля возможных пространственно-временных реализаций гидрометеорологических полей в рассматриваемой области. Необходимый анализ статистической структуры полученного климатического ансамбля показывает, что модель с достаточно высокой точностью воспроизводит корреляционную структуру реальных климатических полей.
Предложенный метод позволяет построить полный набор рассматриваемых гидрометеорологических полей в регулярных пространственно-временных точках, удовлетворяющих климатическим свойствам и соответствующих численной модели динамики атмосферы, т, е, численно восполнить недостаток априорных сведений об этих полях, С другой стороны, предложенный подход с использованием моделирования климатического ансамбля также позволяет с достаточной статистической достоверностью получать оценки свойств результатов численного моделирования.
Рассматриваемая методика иллюстрируется в данной работе на примере исследования климатических выбросов в атмосфере для некоторой локальной области. Сущность этого метода заключается в следующем. Мы предполагаем, что статистическая структура рассматриваемых полей известна приближенно и задана корреляционной матрицей Я, Согласно этой структуре с использованием известных методов статистического моделирования строится ансамбль реализаций этих полей
т
Здесь £(п) = (иг (Х5-, ), Тг (Х^, ), Нг (Х^, 4),...) — вектор реализаций полей скорости, температуры, геопотенциала и т.д. в пространственно-временных точках (Х^, рассматриваемой области; п — размерность этого вектора.
Реализации полей этого ансамбля используются как входные данные для задачи вариационного усвоения с помощью динамической модели, В результате мы получаем новый ансамбль, в котором каждая реализация удовлетворяет динамической модели и статистическая структура полученного нового ансамбля близка к начальной структуре в пределах точности решения задачи вариационного усвоения,
В данной работе рассматривается один из методов статистического моделирования, основанный на спектральном разложении корреляционной матрицы. Для этого представим
Я
Я = Ш ЛШт, (2)
где Ш — матрица собственных векторов корреляционной матрицы Я, а Л — диагональная матрица соответствующих собственных значений. Заметим, что представление (2) являет-
Я
определения главных факторов. Дальнейший шаг в нашем случае состоит в определении
Я
Я1/2 = Ш Л1/2Шт, где
Л1/2
— диагональная матрица, на ЯТ
деляет операцию транспонирования. Тогда метод статистического моделирования можем
определять как
= В€ Я1/2^(г) (Х3 ,у3 Л )Т + ,у3 ^), (3)
где ф(г\хз)т, г = 1, 2,..., — гауссовский случайный вектор с единичной дисперсией и нулевым средним; В^ — диагональная матрица дисперсий; , yj, tj) — вектор оценки среднего для моделируемой стохастической величины. Нетрудно видеть, что матрица корреляций стохастического вектора (3) в точности совпадает с матрицей Я,
Метод вариационного усвоения основан на решении проблемы минимизации функционала качества [3] с использованием метода градиентного спуска. Функционал качества характеризует меру различия между реализацией стохастического поля и поля, полученного решением системы динамических уравнений. Для определения градиента функционала качества используются решения основных и сопряженных задач, соответствующих рассматриваемой модели гидротермодинамики, В конечном результате полученная средняя величина функционала качества, очевидно, определяет меру различия начальной статистической структуры и статистической структуры, полученной при совместном моделировании.
Рассматриваемый метод вариационного усвоения информации с помощью математической модели основан на методе Лагранжа поиска стационарных точек функционала качества при ограничениях, заданных в виде равенств, определении дифференциала Гато, а для построения градиента используется основополагающая идея Г. 11. Марчука, конструктивно предложенная в работе [4] в применении к задачам гидротермодинамики атмосферы и океана.
Заметим, что для эффективного использования вариационного усвоения необходимо теоретическое или численное исследование сходимости и исследование функционала Гато, Кроме того, при построении климатического ансамбля реализаций задачу вариационного усвоения следует рассмотреть на всем интервале моделирования в пределах предсказуемости численной модели, поскольку использование так называемого последовательного усвоения не обеспечивает необходимую гладкость решений и тенденцию для дальнейшего использования полученных полей в режиме прогноза.
Сущность метода вариационного усвоения состоит в следующем. Рассматриваемую численную модель запишем в операторном виде
д Ф
— + А(У,Ф)Ф = 0, (4)
где Ф — вектор состояния; V = Ф^=0 — вектор параметров; Ф)[] — конечно-разностный оператор, определенный системой уравнений рассматриваемого процесса и граничных условий в области С = С ® [0,£], Система уравнений (4) определяет некоторый набор решений, зависящих от вектора параметров V, Проблема состоит в том, чтобы среди этого набора решений выбрать наиболее близкое к конкретной реализации из (3) в смысле некоторого функционального качества, который мы запишем как
к
ГД6 ( , )г>а — скалярное произведение в подпространстве измеренных данных Ф|, в качестве которых используются значения поля размерности в точках по времени ¿к ансамбля реализации (3); Ь — оператор интерполяции решения проблемы (4) в
пространственно-временные точки измерений. Таким образом, необходимо найти минимум функционала .То относительно вектора параметров V при ограничениях (4),
Для решения этой проблемы используется итерационный метод градиентного спуска градиента, основанный на методе Лагранжа и решении прямых и сопряженных задач. Заметим, что в общем случае минимум функционала не единствен, что определяется нелинейностью системы (4) и степенью полноты усвояемых данных, В этом случае минимум функционала определяется начальным приближением итерационного процесса.
Детальное описание объединенной динамико-стохастической модели и ее некоторые характеристики приведены в работе [3],
1. Моделирование климатического ансамбля для локальной области и исследование на его основе климатических выбросов в атмосфере
Таким образом, в соответствии с работами [1, 2] мы рассматриваем климатический ансамбль реализаций многомерных гидрометеорологических полей для выбранного интервала времени и заданной локальной области в виде (1),
Существующие в настоящее время климатические численные модели основаны на полных уравнениях гидротермодинамики, в которых результат моделирования получается за счет все более полного включения физических механизмов, применения детального пространственно-временного разрешения и различных параметризаций (влага, потоки тепла, пограничные слои и т.д.) и интегрирования на достаточно длительном временном интервале с выходом на некоторый квазипериодический режим, В данной работе предлагается моделировать непосредственно статистически независимые реализации пространственно-временных полей для заданной области с оптимально близкой к априорной климатической информации корреляционной структурой, С этой целью используется совместное статистическое и динамическое моделирование, которое позволяет в рамках единого подхода оптимально совместить статистические и гидротермодинамические свойства численного моделирования,
В качестве реальных данных для корреляционной матрицы использованы данные ре-анализа (МСЕР/МСАЕ) полей температуры для зимнего периода за 1948-2005 гг. на десяти стандартных уровнях с дискретностью по времени 6 ч и 2,5x2,5 град, по горизонтали, Выборка осуществлялась для заданной локальной области Северного полушария размером 10x10 град, с центром в точке с координатами 60.56° с.ш. и 77,7° в.д. Задача рассматривается в системе координат ж, у, р в области, нижним основанием которой является прямоугольник в касательной плоскости в центральной точке. Для расчетов выбрано пространственное разрешение 24x20 точек по коорди патам ж и у с шагами Дж = 23.85 км и Ду = 58.74 км,
Я
чтобы с минимальным количеством этих точек, расположенных вне выбранной локальной области, обеспечить соответствующий процесс интерполяции в регулярные узлы расчетной сетки, С учетом десяти уровней по вертикали и двух моментов времени размерность корреляционной матрицы составляет 840,
Заметим, что характер этой матрицы (см, таблицу) качественно соответствует матрице, взятой из работы [5], хотя количественные значения корреляций отличаются, поскольку
Межуровенные коэффициенты корреляций поля температуры, рассчитанные по данным реанализа для выбранной точки локальной области (значения ниже главной диагонали), и климатические коэффициенты [5] (значения выше главной диагонали)
Р, мб Р, мб
1000 850 700 500 400 300 250 200 150 100
1000 1.00 0.67 0.57 0.47 0.45 0.34 - -0.27 - -0.14
850 0.80 1.00 0.74 0.68 0.57 0.29 - -0.45 - -0.45
700 0.70 0.91 1.00 0.76 0.67 0.44 - -0.49 - -0.46
500 0.58 0.78 0.91 1.00 0.94 0.53 - -0.56 - -0.61
400 0.52 0.69 0.82 0.94 1.00 0.67 - -0.55 - -0.70
300 0.23 0.17 0.21 0.32 0.51 1.00 - -0.02 - -0.46
250 -0.13 -0.34 -0.38 -0.34 -0.23 0.63 1.00 - - -
200 -0.21 -0.40 -0.44 -0.43 -0.37 0.35 0.90 1.00 - 0.51
150 -0.17 -0.27 -0.28 -0.25 -0.20 0.35 0.78 0.93 1.00 -
100 -0.09 -0.15 -0.13 -0.10 -0.06 0.36 0.66 0.79 0.93 1.00
в таблице они приведены только для одного из вертикальных профилей температуры в выбранной локальной области,
Я
ций (1). В нашем случае смоделированный ансамбль (1) содержит только пространственно-временные случайные поля температуры, а для получения статистических оценок оказалось достаточным выбрать число реализаций в количестве 2000,
Дальнейший шаг построения климатических реализаций — это решение задачи вариационного усвоения, для которой входными данными являются реализации из (1), В результате мы имеем новый ансамбль реализаций
{&),г =1, 2,...}, (5)
который уже содержит полный набор гидрометеорологических переменных, согласованных относительно численной модели (4) и имеющих корреляционную матрицу для поля
Я
Для анализа динамических процессов реализаций климатического ансамбля (4) будем использовать задачу численного моделирования процесса переноса пассивной примеси в атмосфере, поскольку этот процесс характеризует прежде всего свойства атмосферных движений, и в силу линейности используемой модели переноса примеси характер ее поведения полностью определяется полями скорости ветра климатического ансамбля. Поэтому свойства полученных распределений этой примеси характеризуют также соответствующие свойства климатического ансамбля метеорологических элементов,
В отличие от существующих в настоящее время численных моделей распространения примеси в атмосфере, в которых в качестве фоновых скоростей ветра используются средние климатические поля или поля, выбранные согласно некоторому характерному сценарию исходя, как правило, из субъективных критериев, в настоящей работе предлагается моделировать целый ансамбль реализаций пространственно-временных полей примеси согласно выбранному ансамблю климатических реализаций полей скорости ветра. Это гарантирует в некотором смысле статистическую достоверность получаемых оценок.
Под климатическим распределением примеси в локальной области будем понимать соответствующее среднее по полученному ансамблю реализаций полей примеси.
Для определения понятия климатической траектории используем известный метод индикаторных функций и результат осреднения ансамбля реализаций полей примеси, С этой целью рассмотрим индикаторную функцию в виде
где f — значение тестируемой функции; lev — некоторое пороговое значение. Значение этой функции, равное единице, в рассматриваемой области от усредненной последовательности реализации ансамбля примеси определяет приближенную климатическую траекторию Tr
lev
конкретная траектория Tr» для г-й реализации ансамбля полей примеси,
В численных экспериментах начальное поле примеси, определяющее мгновенный источник загрязнения, выбрано в некоторой области размером 5x5 точек на уровне p = 850 мб с максимальным значением в центре, равным 1, Пороговое значение lev выбрано равным 0,05, Для численного моделирования процесса распространения примеси в работах [1, 2] использовано трехмерное линейное уравнение переноса с аппроксимацией с помощью модификации известной квазимонотонной численной схемы Ван Леера, включенной также в качестве отдельного блока в модель (4), В данной работе использована полулагранжева численная схема [6], При ее реализации использована трехмерная бикубическая сплайн-интерполяция.
Заметим, что оценка степени принадлежности конкретной траектории Tr» распреде-
Tr
в зависимости от целей исследований. Например, оценка принадлежности как отношение числа точек пересечения с климатической траекторией к общему числу точек этой климатической траектории в случае, если конкретная траектория полностью находится в климатической траектории, но остается в достаточно малой окрестности от источника загрязнения, будет иметь малое значение оценки пересечения, В этом случае с точки зрения исследования распространения примеси следует ожидать повышенное содержание примеси в этой области, и мы будем определять данную ситуацию как один из случаев реализации статистического выброса,
С другой стороны, если мы хотим иметь оценку степени геометрической принадлежности конкретной траектории Tr» распределения примеси климатической траектории Tr, необходимо рассмотреть ее в виде отношения числа точек пересечения с климатической траекторией к общему числу точек траектории Tr». Далее будем рассматривать именно эту оценку, обозначив ее через а, Усредненная по времени климатическая траектория, вычисленная по ансамблю полей примеси, представлена на рис, 1, Поскольку реализации полученного ансамбля — стохастические векторные функции, оценка степени принадлежа
стохастической функцией. Вид этой функции приведен на рис, 2, где по оси абсцисс от-
а
0,162, Это указывает на то, что в рассматриваемой задаче процессы распространения примеси имеют некоторое преобладающее направление, и позволяет выделить наиболее вероятную область данного уровня загрязнения. На рис, 2 представлена гистограмма это-
а
приведена осредненная по времени траектория распределения примеси с минимальной степенью принадлежности к климатической траектории, равной 0,193, т.е. ее большая часть находится вне климатической траектории и она может в нашем случае рассматриваться
если f > lev, если f < lev,
(6)
Рис. 1. Осредненные по времени киматическая траектория Тг и траектория распространения примеси, составляющей выброс (5); точка А соответствует центру мгновенного источника примеси; вдоль осей координат отложены номера сеточной области.
Рис. 2. График стохастической функции а и гистограмма ее распределения.
Рис. 3. Оценка вероятности распределения величины а.
B
A
целом результаты, полученные в данной работе, качественно соответствуют результатам
из работ [1, 2], хотя численные оценки несколько различаются. Например, наблюдаются
а
максимум в окрестности этого среднего значения.
Анализ численных результатов показывает, что предложенная техника может использоваться для определения и анализа прямых и обратных климатических траекторий динамических процессов в атмосфере и исследования этих выбросов.
Список литературы
[1] Протасов A.B., Огородников В.А. Численный метод стохастического моделирования климатических траекторий распространения примеси в атмосфере // Тр. Междунар. конф. по вычисл. математике. Ч. I. Новосибирск: Изд-во ИВМиМГ СО РАН, 2004. С. 310-315.
[2] Протасов A.B. Динамико-вероятностное моделирование климатических полей метеоэлементов в локальной области на основе данных реанализа // Тр. Междунар. научно-техн. выставки-конгресса "ГЕО-Сибирь-2006". Т. 3, ч. 2. Новосибирск: Изд-во СГГА, 2006. С. 112— 118.
[3] Ogorodnikov V.A., Protasov АЛ7. Dynamic probabilistic model of atmospheric processes and variational methods of data assimilation // Russ. J. Numer. Anal. Math. Modelling. 1997. Vol. 12, N 5. P. 461-479.
[4] Марчук Г.И. Численное решение задач динамики атмосферы и океана. Л.: Гидрометеоиздат, 1974.
[5] Гандин J1.C., Каган P.J1. Статистические методы интерпретации метеорологических данных. Л.: Гидрометеоиздат, 1976.
[6] Van Loon M. Tests on Semi-Lagrangian Transport and Interpolation. Report NM-R9217 September. CWI is the Centre for Mathematics and Computer Science of the Mathematical Centre Foundation, 1992.
Поступила в редакцию 9 ноября 2006 г.
