CC BY
42
УДК 681.322.05
ПРОГРАММНО-МЕТОДИЧЕСКОЕ ОБЕСПЕЧЕНИЕ ДЛЯ ОБРАБОТКИ И МОДЕЛИРОВАНИЯ АТМОСФЕРНЫХ КЛИМАТИЧЕСКИХ ВЕЛИЧИН
С.Ю.Золотов
В статье рассказывается о новом методе восстановления пропущенных значений климатических величин. Данный метод позволяет строить и прогнозы значений временных рядов соответствующих величин. Также в статье приводится описание модифицированной комбинированной методики и программного обеспечения для расчетов вертикальных потоков тепловой радиации атмосферы.
Инструментальные наблюдения глобальной температуры приземного воздуха, выполненные в течение двадцатого столетия, позволили оценить темп повышения температуры величиной 0,6 °С за сто лет, причем в отдельных регионах (к которым относится и юг Западной Сибири) в последние десятилетия процесс потепления заметно ускорился [1-3]. В качестве основной причины глобального потепления в настоящее время рассматривается увеличение содержания в атмосфере газов, обладающих сильным парниковым эффектом: С02 и СН4.
В зависимости от степени учета соответствующих климатообразующих процессов различные климатические модели дают эффект глобального потепления к 2100 г. от 1,5 до 5,8 °С [1, 2]. Во всех без исключения климатических моделях критически важное значение имеет качество так называемого радиационного блока - совокупности алгоритмов, позволяющих рассчитывать приходящую коротковолновую солнечную радиацию и уходящую от системы Земля-атмос-фера длинноволновую.
Из-за очень жестких требований к быстродействию радиационных блоков климатических моделей в них используются исключительно параметрические и упрощенные методы расчета потоков атмосферной радиации. По мнению ряда специалистов, неудовлетворительная надежность климатических прогнозов обусловлена в значительной мере именно низкой точностью описания радиации в климатических моделях. Таким образом, возникает необходимость в параметризации радиационных процессов с приемлемой точностью.
Для различных вычислений в радиационном блоке климатических моделей необходимо знание численных значений соответствующих климатических величин. Получение экспериментальных данных опирается на ведущиеся на многих станциях мировой сети инструментальные наблюдения за температурой, давлением, влажностью и другими величинами. Инструментальные наблюдения ведутся и за значениями концентрации парниковых газов (углекислый газ, метан, озон) [4]. При этом повышенные требования предъявляются к качеству рядов наблюдений. Анализ имеющихся рядов наблюдений показал, что в них нередки пропуски в силу организационных и технических причин. Таким образом, в непрерывном ряде наблюдений образуются разрывы различной длительности. Для заполнения этих разрывов требуется разработка алгоритмов восстановления пропущенных значений, не снижающих существенно качество всего ряда.
В течение последних двух десятилетий отмечен значительный прогресс в развитии глобальных климатических моделей (ГКМ), обусловленный как достижениями в исследованиях собственно климатической системы, так и ростом возможностей вычислительной техники,
обеспечивающей все большую детализацию и полноту модельных описаний климатически значимых процессов. Современные ГКМ, используемые в расчетах изменений климата, включают в качестве основных компонентов интерактивные (взаимодействующие друг с другом) модели атмосферы, океана, верхних слоев суши, криосферы и биосферы. Однако их применение в расчетах изменения климата является довольно сложной задачей, и поэтому в настоящее время все еще актуальны прогнозы, основанные на использовании рядов непосредственно самой климатической величины.
Восстановление/прогноз значений климатических величин
с помощью вейвлет-преобразования
Ряд климатических величин измеряется в дискретные моменты времени и в силу различных причин имеет пропуски. Обычно для восстановления пропущенных значений применяется определенная математическая аппроксимация временного ряда (различные полиномы, сплайны). Специфика всех этих аппроксимаций такова, что на практике ищут коэффициенты полинома/сплайна по очень малому количеству близлежащих точек и, таким образом, не учитывается информация о колебательной структуре всего ряда.
В исследованиях колебательной структуры прежде всего используют преобразование Фурье. Преобразование Фурье является полезным математическим аппаратом для анализа сигналов, однако иногда оказывается недостаточно эффективным при обработке сложных сигналов. Оно не позволяет анализировать локальные свойства сигнала, так как базисные функции Фурье-преобразования определены на всей временной оси. Поэтому в последнее время все чаще вместо Фурье-преобразования используется вейвлет-преобразование.
Вейвлет-преобразование сигнала состоит в его разложении по базису, построенному из обладающей определенными свойствами специальной функции (вейвлета) посредством масштабных изменений и переносов. В отличие от преобразования Фурье вейвлет-преобразование обеспечивает двухмерный анализ исследуемого одномерного сигнала, при этом частота и время рассматриваются как независимые переменные.
На рис. 1 приведены вейвлет-коэффициенты для модельного ряда ) , tе[0, 0,1п, ...
..., 18п]. На этом рисунке видна четкая масштабная повторяющаяся структура, что дает нам основание с достаточной точностью интерполировать в местах пропуска, а также продолжить по временной оси эти вейвлет-коэффициенты. Далее с помощью обратного вейвлет-преоб-разования восстанавливается ряд. Результаты восстановления в точках пропуска ^ е [10п, 10,1п, ...,12п] функции вт(/) , заданной на интервале tе[0, 0,1п, ..., 18п], и прогноза на моменты времени I е [18,1п, 18,2п,..., 20п] сведены на рис. 2.
Для сравнения различных методов интерполяции возьмем ряд среднегодовых значений температуры г. Омска в период 1916-2002 гг. Интерполяция проведена на период 1956-1959 гг. Из всех известных методов интерполяции для сравнения были выбраны методы интерполяции полиномом Ньютона третьей степени и кубическим сплайном как широко распространенные и встроенные во все популярные математические пакеты. Анализ результатов интерполяции стандартными методами показывает, что колебательная структура ряда полиномом/сплайном абсолютно не учитывается. Как следствие, наблюдается большое отклонение восстановленных стандартными методами данных от реальных.
Проведем интерполяцию с помощью вейвлет-преобразования. Анализ карты абсолютных значений вейвлет-коэффициентов показывает четко выраженные колебания масштабов 11 и 28 лет, периодически проявляется масштаб в 5,5 лет (рис. 3).
2я 4л бя 8л 10* 12 я 14л ]бя 18л Время
Рис. 1 - Карта абсолютных значений вейвлет-коэффициентов ряда 5т(Т)
2л 4л 6л 8л Юл 12л 14л 16л 18л 20л
Бремя
Рис. 2 - Восстановление ряда $т(Т) (пунктир) в диапазоне t е [10П, 10,1П,..., 12П] и прогноз (пунктир) на / е [18,1п, 18,2п, ..., 20п]
1920 1930 1940 1950 1950 1970 1930 1990
Время, г.
Рис. 3 - Карта абсолютных значений вейвлет-коэффициентов годового температурного ряда г. Омска
Далее мы восстанавливаем пропущенные вейвлет-коэффициенты и по ним рассчитываем значение температурного ряда в точках пропуска данных, т.е. получается непосредственно интерполяция.
Сравнительный анализ различных методов интерполяции приведен на рис. 4. Из него видно, что метод интерполяции на основе вейвлет-преобразования имеет преимущество по сравнению с другими интерполяционными методами.
2_I_I_I_I_I_I_I_I_I_I_I_I_I_
1950 1952 1954 1956 1958 1960 1962 1964
Время, г.
Рис. 4 - Сравнение различных методов интерполяции (вейвлет-преобразование - пунктир, полином Ньютона -точки, кубический сплайн - штрих-пунктир) с реальными данными (сплошная линия) среднегодовой температуры г. Омска (1916-2002 гг.). Вертикальными отрезками показаны 90 %-е доверительные интервалы
Моделирование атмосферных тепловых потоков
Ключевую роль в изменении климата играют радиационные процессы, так как они прямо или косвенно связаны с такими климатообразующими факторами, как вариации солнечной постоянной, трансформации свойств подстилающей поверхности, изменения газового и аэрозольного состава атмосферы и др. Одной из проблем моделирования радиационных процессов является учет молекулярного поглощения газовыми составляющими атмосферы вследствие высокой селективности коэффициентов поглощения.
Для расчета радиационных потоков необходимо знать распределение метеорологических величин вдоль трассы луча. Обычно при проведении расчетов потоков излучения атмосферу считают послойно однородной. Такое допущение во многих случаях вполне оправдано и, что особенно важно, позволяет уменьшить объем метеорологической информации путем замены распределения метеопараметров вдоль произвольных направлений на распределение лишь вдоль вертикального направления с последующим перерасчетом на заданный угол трассы.
В связи с тем что реальные измерения высотного профиля некоторой метеорологической величины для любого географического района и в любое время провести практически невозможно, на практике используется следующая информация.
1. Высотные распределения температуры, давления и концентрации газов в виде статистических моделей, адекватно описывающих данный географический район и сезон года. В настоящее время активно используются полученные в Институте оптики атмосферы СО
РАН (г. Томск) статистические модели высотного распределения термодинамических параметров. Эти модели содержат не только средние значения температуры, давления и концентрации газов, но и параметры их изменчивости в пространстве и во времени (среднеквадратические отклонения и матрицы ковариации).
2. Высотные распределения давления, температуры и концентрации газов, полученные косвенным путем из измерений оптических характеристик атмосферы различными приборами (лидарами, спектрофотометрами, озонометрами).
3. Реальные значения метеовеличин у поверхности земли и восстановленные на их основе с помощью статистических методов высотные профили.
Таким образом, рассчитанный профиль получается более близким к реальному по сравнению с модельным (см. пункт 1) за счет использования реальных значений метеовеличин у поверхности земли. Для расчетов необходимо знать основные статистические значения величины по высоте (среднее, среднеквадратическое отклонение, корреляционную матрицу). Обычно рассматриваются три типичных для данной ситуации случая: а) известно значение величины q на высоте к0 = 0; б) известно значение искомой величины q и другой, коррелированной с
ней, величины и на высоте к = И00; в) известно значение величины q на высотах к0 = 0 и
Их > к0.
В качестве примера возьмем высотные профили (0-30 км) измеренных значений температуры и удельной влажности (0-8 км) над г. Новосибирском в период с 1 июля 1961 г. по 24 августа 1970 г. (483 временных точки).
На рис. 5 приведено сравнение реального температурного профиля за 26 июля 1965 г. и оценки профиля, рассчитанного по вышеописанному случаю (б). Из рис. 5 видно, что достигается хорошее соответствие между оценкой и реальным результатом.
Рис. 5 - Сравнение реального (сплошная линия) температурного профиля за 26 июля 1965 г., модельного (точки) (лето средних широт) и расчетного (пунктир) для случая (б) (см. текст). Горизонтальными черточками показана расчетная ошибка восстановления профиля
Для моделирования распространения теплового излучения Земли используют уравнение переноса радиации. Основной характеристикой этого уравнения является функция пропуска-
ния. Именно для ее расчета требуются большие затраты машинного времени. В связи с этим многие исследователи пытаются создать эффективные модели расчета характеристик молекулярного поглощения, которые позволяли бы обеспечить высокую степень параметризации без потери точности.
В зависимости от решаемых задач оптики атмосферы функции пропускания должны быть известны с высоким (например, для интерпретации данных спутниковых приборов) либо грубым (например, для решения радиационно-климатических задач) спектральным разрешением. В настоящее время для расчета функций пропускания со средним и низким спектральным разрешением широко используются следующие группы методов: методы на основе полинейного счета (Ипв-Бу-Ипв: ЬБЦ метод ^-распределения); эмпирические методы (методика ГОИ, методика ГИПО, методика Ю^РАЫ, методика ГГО); методы на основе моделТаскрегтоа екш© щедеяи полос поглощения основаны на идеализированном описании спектров поглощения. Эти модели - параметрические, причем сами параметры определяются с помощью экспериментальных данных или рассчитываются прямым методом для различных значений термодинамических параметров среды. На основе моделей полос поглощения была разработана комбинированная методика расчета функций пропускания. По этой методике каждая модель применяется к тому спектральному интервалу, в котором она обеспечивает наилучшую точность. Достоинствами комбинированной методики являются: а) быстрый расчет пропускания (параметры для каждой длины волны считаются постоянными величинами); б) возможность проводить аналитические вычисления при решении радиационных задач. К недостаткам комбинированной методики следует отнести отсутствие границ применимости и небольшую точность. Необходимо отметить еще и такой факт, что параметры в разных диапазонах длин волн найдены при разных спектральных разрешениях.
Для устранения больших погрешностей при расчете функции пропускания нами была создана модифицированная комбинированная методика (МКМ). Ее преимущества:
а) подбор всех параметров модели проводится с одним спектральным разрешением;
б) расчет параметров для всех газов осуществляется в широком диапазоне волновых чисел (в практической реализации 330-10000 см-1);
в) варьирование модельных параметров производится в зависимости от значения поглощающей массы газа.
Определение параметров некоторой модели осуществляется с помощью метода наименьших квадратов. Исходные данные получаются либо из эксперимента, либо при помощи метода ЬБЬ. Для дальнейшего улучшения точности расчетов максимально возможный диапазон поглощающих масс разбивается на несколько поддиапазонов.
В таблице приведено сравнение точности расчета пропускания МКМ и Ю^РАЫ-7 с ЬБЬ. Для расчетов были выбраны следующие условия: спектральное разрешение 20 см-1, вертикальная трасса 0-100 км, зенитный угол равен нулю, спектральный диапазон 330-10000 см-1 с шагом 5 см-1 . Для подбора параметров моделей в модифицированной комбинированной ме-
тодике общий диапазон поглощающих масс был разбит на десять поддиапазонов. Из таблицы видно, что МКМ гораздо точнее, чем методика Ю^РАЫ.
Сравнение расчетов функции пропускания по МКМ и LOWTRAN-7 с 1_В1_ (лето средних широт)
Газ Макс. погреш. МКМ, % Сред. погреш. МКМ, % Макс. погреш. LOWTRAN-7, % Сред. погреш. LOWTRAN-7, %
H2O 10,1 0,6 54,0 7,9
CO2 34,0 0,5 86,2 1,7
O3 0,2 меньше 0,1 29,1 0,3
n2o 3,3 меньше 0,1 20,4 0,5
CO 0,9 меньше 0,1 10,4 0,1
CH4 6,1 0,1 30,5 1,1
O2 0,4 меньше 0,1 35,1 0,7
NO меньше 0,1 меньше 0,1 0,1 меньше 0,1
SO2 2,2 меньше 0,1 68,1 0,5
no2 4,0 меньше 0,1 3,9 меньше 0,1
NH3 0,1 меньше 0,1 4,2 0,1
HNO3 меньше 0,1 меньше 0,1 - -
Программное обеспечение
Для расчетов вертикальных потоков тепловой радиации атмосферы с помощью МКМ разработан программный комплекс «ABS» (ABSorption). В его состав входят: а) ПО расчета потоков эмиссии метана; б) ПО прогноза и восстановления пропущенных значений климатических величин с помощью вейвлет-преобразования; в) ПО модифицированной комбинированной методики (рис. 6). Особенность данного программного комплекса заключается в возможности отдельного использования каждого ПО, входящего в его структуру.
Рис. 6 - Структура программного комплекса «ABS»
Программное обеспечение было создано для реализации оригинального авторского алгоритма прогноза и восстановления пропущенных значений климатических величин с помощью вейвлет-преобразования. Программное обеспечение реализовано в системе MATLAB.
Алгоритм прогноза и восстановления пропущенных значений климатических величин с помощью вейвлет-преобразования следующий:
1. Считывание входного временного ряда.
2. Считывание вида действия, которое необходимо осуществить: восстановить пропущенные значения во времени или сделать прогноз. В случае выполнения процедуры интерполяции необходимо ввести временные интервалы для восстановления, а если выбран прогноз -то количество временных точек для данного прогноза.
3. Расчет коэффициентов вейвлет-преобразования.
4. Восстановление вейвлет-коэффициентов во время пропуска или их продолжение по временной оси.
5. Восстановление исходного ряда из измененных вейвлет-коэффициентов. Автоматически в ходе данной процедуры получается интерполирование временного ряда или его прогноз.
6. Вывод в виде графика (эта возможность предусмотрена опциально):
а) входной временной последовательности;
б) результата работы программы (т.е. восстановленных значений или прогноза);
в) карты абсолютных значений вейвлет-коэффициентов.
7. Возвращение пользователю:
а) в виде выходного ряда - полученный результат;
б) в виде матрицы - полученные вейвлет-коэффициенты.
Программное обеспечение обработки данных эмиссии метана для установки, описанной в работе [4], до этого не существовало, так как вся обработка проводилась в виде определенных запросов на SQL, реализованных в среде MS Access. Авторское программное обеспечение для расчета потоков метана было создано в системе MATLAB.
Алгоритм поиска значений потоков эмиссии метана следующий:
1. Считывание данных из входного файла датчика метана и их запись в структуру входных данных. По умолчанию каждому входному файлу соответствует собственная структура, однако при необходимости эти структуры могут быть объединены. Основу входной структуры составляют записанные во времени значения напряжения метанового детектора. Данная структура является основным входным источником данных для авторского программного обеспечения.
2. Использование априорной информации об искомых величинах, например, задание условия, что величина эмиссии потока метана не будет превышать некоторого значения в данный период времени.
3. Игнорирование временных отрезков, в которых самим микропроцессором зафиксированы сбои оборудования. Таким образом, поиск неизвестных величин осуществляется только по корректно регистрируемым значениям напряжения детектора CH4 в текущем часе.
4. Выбор пользователем временных подынтервалов, в которых необходимо осуществить подбор неизвестных величин. Эта процедура позволяет игнорировать заведомо известные временные отрезки, в которых особенно проявляются как нестабильность работы аппаратуры, так и сильные помехи. В отличие от пункта 3, данные временные подынтервалы не считаются микропроцессором как сбои оборудования, и только сам пользователь должен решать, как сильно они будут искажать реальные значения искомых величин. Для этого в помощь пользователю строится специальный график временного хода значения напряжения детектора метана, по которому он и принимает соответствующее решение. Данная процедура является опци-альной и может быть отключена (т.е. поиск неизвестных величин будет проходить по всем зарегистрированным значения напряжения детектора).
5. Поиск неизвестных величин с использованием метода наименьших квадратов. Для этого используется встроенная в MATLAB специальная функция поиска минимума целевой функции с ограничениями.
6. Строится график как реальных значений напряжения, так и расчетных. Такой вывод результата позволяет косвенно судить о том, насколько близки к истинным подобранные в пункте 5 искомые величины. Данная операция носит опциальный характер.
7. Результаты поиска неизвестных величин записываются в выходную структуру. Основу данной выходной структуры составляют записанные во времени (за каждый час) значения подобранных величин. При необходимости эти величины могут быть показаны на графике.
Программное обеспечение МКМ представляет собой динамически подключаемую библиотеку и позволяет рассчитывать пропускание/поглощение, восходящую или нисходящую интенсивность теплового излучения, потоки тепловой радиации атмосферы.
ЛИТЕРАТУРА
1. Climate Change 2001. The Scientific Basis. Contribution of Working Group I to the Third Assessment Report of the IPCC. Summary for Policymakers and Technical Summary. - WMO/UNEP, 2001.
2. Израэль Ю.А. Изменения глобального климата. Роль антропогенных воздействий / Ю.А. Израэль, Г.В. Груза, В.М. Катцов, В.П. Мелешко // Метеорология и гидрология. - 2001. -№ 5. - С. 5-21.
3. Региональный мониторинг атмосферы. Ч. 4 / Под ред. д-ра физ.-мат. наук чл.-кор. РАН М.В. Кабанова. - Томск: МГП «Раско», 2000. - 270 с.
4. Золотов С. Ю. Потоки метана заболоченной территории Западной Сибири / С.Ю. Зо-лотов, А.И. Надеев // Доклады Томского государственного университета систем управления и радиоэлектроники. Т. 8. Автоматизированные системы обработки информации, управления и
проектирования: Сб. науч. трудов. - Томск: Том. гос. ун-т систем упр. и радиоэлектроники, 2003. - С. 17-20.
