Материал: Моделирование колебаний климата

В работе [1] В.И. Бышев и И.В. Серых с соавторами выдвигают гипотезу о существовании моды глобальной атмосферно-океанической циркуляции, охватывающей весь земной шар и передающей влияние от Эль-Ниньо в Тихом океане в Атлантику, Индийский океан, Северный-Ледовитый океан и по всему земному шару. Анимации представленные в работе [26] не противоречат возможности распространения глобальных мод климатических колебаний, в частности Эль-Ниньо, по всей атмосфере Земли. Можно предположить, что ввиду сильного перемешивания вод северной Атлантики и наличия в ней мощных течений, таких как Гольфстрим, ее влияние на климат Земли в целом довольно велико. По результатам предварительного анализа Международных климатических моделей CMIP5, проведенного И.В. Серых, они плохо воспроизводят МАК. Тем не менее, при надежных начальных условиях, МАК воспроизводится и охватывает океан до больших глубин [33].

2.2 Модель авторегрессии

Модель авторегрессии (АR) задается формулой

где  - параметры модели,  - поступающий на вход белый шум,- порядок модели[5].

Пусть имеются экспериментально полученные уровней временного ряда, заданные соотношением

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

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

2.3 Метод средней квадратической коллокации

,

Теперь введем матрицы (ковариационные)

где  и  - автоковариационные матрицы (АКМ) векторов и соответственно, - АКМ между и .Тогда требуется оценить вектор -сигнала, при известных АКФ и векторе данных . После этого осталось искать оценку  -сигнала следующим образом

где  -матрица . Или же каждое значение  - вектора аппроксимируется с вектором - данных.

АКМ вектора ошибок (определяемая соотношением )

Диагональные элементы АКМ являются дисперсиями значений сигнала

Задача фильтрации временных рядов:

Пусть дана последовательность оценок некоторого стохастического  - сигналав виде -мерного вектора. Тогда,  - сигнал в виде -мерного вектора , смешанного со случайными коррелированными ошибками , т.е

Если известна априорная АКМ шума , то задача фильтрации/выделения сигнала, при :

где -АКМ данных . Для вычисления ковариационной матрицы  достаточно найти АКФ ряда . В этом случае  Вычисляем матрицу

Далее осталось лишь найти АКМ сигнала

Если же иметь дело с белым шумом то модель принимает иной вид

Метод СКК позволяет осуществлять прогнозирование временных рядов с помощью АКФ. Смещенная оценка которой вычисляется так

Значения АКФ заносятся в симметрическую ковариационную матрицу , так что .

Задача прогноза временных рядов:

Прогноз, в свою очередь, вычисляется как

,

где -часть матрицы  , а -часть матрицы[3].

Для заполнения ковариационной матрицы эмпирическими оценками АКФ в книге [4] предложено предварительно отфильтровать значения АКФ, что позволяет гарантированно симметризовать АКМ и обеспечить ее положительную определенность.

Алгоритм фильтрации оценки ковариационной функции  таков:

(а) путем локальной аппроксимации АКФ  в окрестностях точки  определяем скачок этой функции в нуле , и затем оцениваем спектральную плотность  белого шума;

(б) вычисляем оценку функции спектральной плотности

(в) осуществляем фильтрацию функции ; для этого все значения полагаем равными ; осуществляя другое преобразование - сглаживание, в результате получаем положительную для всех частот функцию спектральной плотности

(г) с помощью обратного преобразования Фурье

находим заведомо положительную оценку АКФ [3].

3.Результаты исследований

.1 Анализ временных рядов колебаний климата

индекс климатический океан

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

Индекс САК содержит разность давлений на станциях Гибралтар и Рейкьявик с 1900 по 2015 годы, с шагом в один месяц.

Исходный ряд САК центра Hurrell [28]сравнивался со сглаженным рядом разности давлений в Исландском и Азорском центрах действия атмосферы по данным [12] (далее ряд Vil’fand). Результаты сравнения представлены на рисунке 8.

Рис. 8. Индексы САК по данным Hurrel (синим) и по данным Vil’fand (красным)

Далее сглаженный ряд САК с годовым шагом был проинтегрирован и поделен на 40 (для соответствия шкал). Результат представлен на рисунке 9в сравнении с индексом МАК.Вычисленная интегральная кривая САК, хорошо согласуется с МАК. Очевидно, что гипотезы Гулева и Пенланд о том, что Атлантика интегрирует САК и переводит его в МАК, имеют под собой основания.

Рис. 9. График МАК и проинтегрированный САК

Индекс МАК с шагом в один год отражает значения аномалий температуры поверхности Атлантического океана севернее экватора с 1850 по 2014 год. На рис. 9, 10мы использовали ряд МАК центра AOML/NOAA [29].

Теперь сопоставим индекс МАК с аномалиями температуры по всей Земле (океан+суша) HadCRUT4 (см. рис. 5 слева), после снятия параболического тренда.

Из рисунка 10видно прекрасное согласие по фазе и по амплитуде индекса МАК и глобальных аномалий температуры на Земле. Отличие составляет лишь высокочастотная часть аномалий. В целом же, Северная Атлантика оказывает сильнейшее влияние на климат Земли.

Рис. 10. Аномалии глобальной температуре на Земле и индекс МАК

Далее для задач прогнозирования САК выполнялось моделирование его детерменированной и стохастической частей. На первом этапе выполнялось моделирование трендов САК. Полиномиальный тренд второго порядка был подобран с помощью функции polyfit. Гармоника с периодом 55 лет подбиралась метом наименьших квадратов. Модели трендов представлены на рисунке 11. Остаточные разности без трендов моделировались как процессы авторегрессии 10 порядка. Параметры авторегрессии, дисперсия исходного белого шума оценивалась функцией ARMatLab методом Берга. Прогноз САК строился на основе суммирования прогнозов полиномиальной, гармонической и авторегрессионной частей. На рисунке 12 представлены исходные ряды (синим) и их прогнозы (красным) для рядов САК Hurrell (слева)и Vil’fand(справа).

Рис. 11.Модели гармонического (слева) и полиномиального (справа) трендов не сглаженного Hurrell (вверху) и сглаженного Vil’fand (внизу) индекса САК

Рис. 12. Прогностическая модель САК до 2026 года по рядам Hurrell (слева) и Vil’fand (справа)

По той же методике выполнялось моделирование и прогнозирование ЭНЮК. Использован индекс ЭНЮК, построенный поразностямзначений приземного давления между станциями Таити и Дарвин с 1890 по 2015 годы, с шагом в один месяц. На первом этапе для исходного ряда ЭНЮК центра Australian Bureau of Meteorology[30] была построена АКФ и периодограмма. В нашей работе наряду с классическим спектральным анализом [4] используется построение фильтрованной периодограммы по методу Блэкмана-Тьюки [31]. Он базируется на теореме о свертке, согласно которой сглаживание спектра в области частот соответствует умножению АКФ на временное окно, центрированное в нуле.

К примеру, на рисунке 13представлена автоковариационная функция для ЭНЮК. Она имеет значительный скачок вблизи нуля и максимумы, соответствующие квазипериодам 28, 36, 60 месяцев и др.

Рис. 13. Автоковариационная функция ЭНЮК

С использованием Фурье-преобразования АКФ, построенной по исходным данным индекса ЭНЮК, вычислена периодограмма на рисунке 14.

Рис. 14. Спектр ЭНЮК, построенный по исходным данным (синим) и их АКФ (красным)

Сглаженная периодограмма индекса ЭНЮКпо методу Блэкмена-Тьюкис корреляционным окном длиной 33 года представлена на рисунке 15. Четко заметны периоды, кратные периоду Чандлеровского колебания (14 месяцев) и др.

Рис.15. Периодограмма с выделенными пиками на некоторых из периодов

Рис. 16. Гармонический (слева) и полиномиальный (справа) тренды ЭНЮК

Полиномиальный тренд второго порядка был подобран с помощью функции polyfit. Гармоника с периодом 85 лет подбиралась метом наименьших квадратов. Модели трендов представлены на рисунке 16.