Известия РАН. Физика атмосферы и океана, 2023, T. 59, № 6, стр. 686-706

1) Метод k-means

Метод k-means (KM) в настоящее время является одним из самых распространенных и простых методов кластерного анализа, который позволяет разделять объекты кластеризации из набора данных на хорошо различимые группы (кластеры), выполняя задачу максимизации межкластерного и минимизации внутрикластерного расстояния лучше большинства других методов.

В случае выделения погодных режимов, объектами кластеризации являются суточные (или среднемесячные) поля атмосферной циркуляции, которые можно представить как точки, заданные в пространстве размерностью L. L будет определяться либо числом узлов сетки, на котором заданы поля, либо числом первых главных компонент разложения полей на эмпирические ортогональные функции, на которые предварительно, т.е. перед процессом кластеризации, были разложены поля, что, благодаря снижению размерности, позволяет получить более устойчивые результаты. Далее в тексте, представленные таким образом суточные поля циркуляции будут называться просто “точки”.

Перед началом работы алгоритма задается число K будущих кластеров и функция расстояния в пространстве набора данных, в качестве которой обычно используется Евклидово расстояние:

где X_i – точка из набора данных (в задаче выделения погодных режимов – оригинальное или разложенное на эмпирические ортогональные функции суточное поле давления и/или высоты геопотенциала), N – число точек в наборе данных, L – размерность пространства, в котором заданы точки.

Алгоритм KM состоит из следующих шагов:

1) Некоторым случайным образом задаются K начальных положений центроидов будущих кластеров.

2) Точки набора данных присваиваются тем кластерам, к центроидам которых они ближе всего.

3) Положения центроидов кластеров пересчитываются как центры точек, присвоенных данному кластеру на предыдущем пункте.

4) Пункты 2 и 3 повторяются до тех пор, пока алгоритм не сойдется, то есть положения центроидов не перестанут меняться, а внутрикластерная дисперсия не перестанет убывать:

где V_iter – суммарная внутрикластерная дисперсия, C_k,iter – центроид k-ого кластера, N_k – число точек, попавших в k-ый кластер, iter – номер итерации алгоритма, itermax – итерация, при которой алгоритм сойдется.

Итоговое положение центроидов зависит исключительно от их начального положения, что, при их неоптимальной инициализации, приводит к возможному схождению алгоритма к локальному оптимуму, для которого итоговое V будет больше, чем могло быть. В связи с этим рекомендуется множественный прогон алгоритма по одному и тому же набору данных и последующий отбор решения с минимальным V [Bradley et al., 1998]. Для улучшения инициализации алгоритма KM существует оптимальный метод задания начальных положений центроидов – “k-means++” – при котором центроиды последовательно задаются с вероятностью, пропорциональной расстоянию до ближайшего из уже заданных центроидов. В работе [Arthur et al., 2007] показано, что такая инициализация центроидов позволяет алгоритму сходится с меньшими V за меньшее число итераций. Алгоритм инициализации “k-means++” состоит из следующих шагов:

1. Первый центроид C₁ выбирается полностью случайно с координатами одной из точек кластеризуемого набора данных.

2. Второй центроид выбирается по координатам оставшихся точек с вероятностью, пропорциональной d²(X_i,C₁) (формула 1). Вероятности нормируются так, чтобы в сумме по всем точкам вероятность равнялась 1.

3. Если K > 2, то все последующие k центроиды (k <= K) выбираются случайно по координатам точек с вероятностью, пропорциональной $\mathop {\min }\limits_{j = 1:k - 1} {{d}_{2}}({{X}_{i}},{{C}_{j}})$, т.е. пропорциональной расстоянию от точки до ближайшего из уже заданных центроидов.

Для устранения детерминированности метода KM, то есть зависимости итогового решения исключительно от инициализации центроидов, существует оптимизация “simulated annealing” (имитация отжига) [Selim et al., 1991]. Во время работы алгоритма на шаге присвоения точек кластерам точки могут быть отнесены не к кластеру ближайшего центроида, а к любому другому с вероятностью P, обратно пропорциональной удаленности от точки до центроида этого кластера:

где $d_{{{\text{old}}}}^{2}$ – расстояние от точки до центроида, к которому она была приписана на предыдущей итерации алгоритма, $d_{{k,{\text{new}}}}^{2}$ – расстояние от точки до k-го центроида, к которому она потенциально будет приписана по итогу текущей итерации алгоритма. T_iter – параметр, убывающий с каждой последующей итерацией – T_{iter + 1} = CT_iter, где C – так называемый коэффициент охлаждения, очень близкий к единице и подбираемый вручную.

Коэффициент охлаждения C подбирается таким образом, чтобы в начале работы алгоритма точки часто приписывались к “неправильным”, т.е. не ближайшим, кластерам с целью обойти решения с локальными оптимумами V, но с ростом числа итераций сходились к итоговому решению, когда вероятности P будут настолько малы, что вероятности быть присвоенной к не ближайшему кластеру перестанут реализовываться, а алгоритм остановится по тому же критерию, что и оригинальный KM (формула 3). Стоит отметить, что при $d_{{{\text{old}}}}^{2}$ = $d_{{k,{\text{new}}}}^{2}$ вероятность P максимальна, то есть вероятность для точки остаться в ближайшем кластере выше (по формуле 4 она равна 1, но эти вероятности нормируются так, чтобы для каждой точки сумма вероятностей быть приписанной к любому из k кластеров равнялась 1), чем в более удаленном.

В работе по кластеризации суточных полей давления в Евро-Атлантике показано [Philipp et al., 2007], что использование оптимизации “simulated annealing” повышает воспроизводимость и лучше минимизирует V итоговых решений по сравнению с традиционным алгоритмом k-means.

2) Иерархическая кластеризация Уорда

Иерархическая кластеризация представляет собой подход, при котором данные разбиваются на систему вложенных кластеров, создавая “дерево” кластеров, которое можно визуализировать с помощью так называемой дендрограммы. Иерархическая кластеризация принципиально разделяется на два подвида – дивизивную (“сверху-вниз”) и агломеративную (“снизу-вверх”) [Roux, 2018].

При дивизивной кластеризации набор данных на начальном этап представляет собой единый кластер и последовательно разбивается на подкластеры, однако такой подход используется редко – для него не существует устоявшихся методов, например, на каждом этапе разбиения можно использовать k-means или другой метод [Lamrous et al., 2006], при этом на каждом этапе исследователю придется каким-то образом выбирать число кластеров K и заново прогонять алгоритм. Агломеративная кластеризация используется чаще [Murtagh et al., 2012], при этом подходе, наоборот, на начальном этапе каждая точка представляет собой отдельный кластер, далее кластеры попарно объединяются по одному из выбранных методов сцепки кластеров (linkage method), которых существует несколько [Murtagh et al., 2017].

На основе выбранного метода сцепки на каждой итерации ищется и объединяется такая пара кластеров, для которой будет минимален параметр, зависящий от выбранного метода сцепки [Govender et al., 2020]: 1) Single linkage, метод ближайшего соседа – расстояние двух ближайших точек из разных кластеров; 2) Complete linkage (метод дальнего соседа) – расстояние двух наиболее удаленных точек из разных кластеров; 3) Average linkage (метод средней связи) – среднее попарное расстояние между точками из разных кластеров; 4) Centroid linkage (центроидный метод) – расстояние между центроидами кластеров; 5) Ward’s linkage (метод Уорда) – расстояние считается как потенциальная прибавка к внутрикластерной дисперсии пары кластеров после их объединения, то есть ищется и объединяется такая пара кластеров (например, A и B), для которой ΔV будет минимально:

где P и Q – любая пара кластеров из имеющихся на данной итерации алгоритма. После объединения, A и B уже не рассматриваются как отдельные кластеры, а как новый кластер AB, алгоритм поиска пары с наименьшим ΔV начинается заново и продолжается до тех пор, пока все кластеры не сольются в один.

Применительно к задаче идентификации погодных режимов, при использовании иерархической кластеризации обычно применяется метод сцепки Уорда [Cheng et al., 1993], который похож на k-means тем, что основным минимизируемым параметром в процессе работы алгоритма здесь является внутрикластерная дисперсия. Преимуществом иерархической кластеризации Уорда (HW) является отсутствие зависимости от инициализации кластеров, которая в принципе отсутствует в HW. Конечный результат кластеризации методом HW зависит исключительно от кластеризуемых данных, при этом по итогу кластеризации исследователь получает возможность изучить дендрограмму кластеризации и остановиться на итерации с интересующим его числом кластеров K. В работе [Bao et al., 2015] кластеризация HW используется для выделения зимних режимов атмосферной циркуляции в Северном полушарии и сравнивается с результатами по методу SOM.

3) Модель гауссовой смеси

Одним из подходов при кластеризации многомерных данных является исследование их плотности распределения. Например, кластеры можно определять путем поиска локальных функций плотности распределения данных. При таком подходе в пространстве кластеризуемого набора данных ищутся локальные максимумы плотности распределения, представляющие из себя точки, имеющие максимальное число соседей (других точек) в некотором радиусе R многомерной сферы по сравнению с другими точками, который подбирается исследователем эмпирически. Такой метод, например, ранее использовался при выделении погодных режимов в Северном полушарии [Molteni et al., 1990]. Применительно к задачам выделения погодных режимов данный метод кластеризации последнее время применяется не так часто, этот метод похож на один из известных методов кластерного анализа под названием DBSCAN (Density-based spatial clustering of applications with noise) [Khan et al., 2014], преимущество которого состоит в способности отделять шум из набора данных, а также малом числе входных параметров перед кластеризацией, всего 2 параметра – радиус R и минимальное число соседей для формирования кластера. Одно из свойств данного метода, которое можно отметить как недостаток – формирование одного крупного, перенаселенного кластера, и большого числа небольших.

Частым предположением о природе кластеризуемых данных служит гипотеза, что данные распределены нормально. В случае кластеризации, когда исследователь полагает, что набор данных состоит из нескольких хорошо различимых групп, можно также предположить, что объекты внутри каждой из этих групп также распределены нормально [Banfield et al., 1993]. Кластеризация гауссовой смесью (GM), в отличие от KM и HW, относится к типу “мягких кластеризаций” (soft clustering) [Kearns et al., 1998], когда по результатам работы алгоритма каждый объект приписывается не к одному определенному кластеру, а к каждому из них с определенной вероятностью, однако исследователю не запрещается по результатам кластеризации присвоить каждый объект из набора данных исключительно тому кластеру, к которому он принадлежит с наибольшей вероятностью.

При кластеризации GM в пространстве набора данных инициализируется K многомерных нормальных распределений, то есть кластеров, для каждого из которых каким-либо способом задаются (угадываются) параметры – среднее и матрица ковариаций [Bilmes et al., 1998]. Среднее при этом можно считать аналогом центроида кластера. Также для каждого распределения задаются веса, которые можно трактовать как доля точек выборки, принадлежащих данному кластеру. Далее, в процессе так называемого EM-алгоритма (Expectation-maximization algorithm) [Yang et al., 2012], основанного на теореме Байеса, и который, вообще говоря, можно использовать и для смеси других распределений, средние, матрицы ковариаций и веса кластеров итеративно пересчитываются. Алгоритм можно представить в виде следующих этапов:

1. Инициализируются средние μ_k, матрицы ковариаций Σ_k и веса кластеров w_k. Значения задаются, как правило, случайным образом, чтобы можно было прогнать алгоритм несколько раз и сравнить результаты. Для инициализации можно использовать, например, способ “k-means++” или найти K первых локальных максимумов плотности распределения точек набора данных по методу, кратко упомянутому в начале данного раздела и более подробнее описанного в [Molteni et al., 1990].

2. Для каждой точки X_i, i = (1…N) из набора данных пересчитывается вероятность ее принадлежности к одному из кластеров (многомерных нормальных распределений):

где B_k – одно из K многомерных нормальных распределений, ${{\rho }_{{{{B}_{k}}}}}$ – функция плотности вероятности этого распределения, для многомерного нормального распределения определяемая по формуле: где L – размерность данных, μ_k – среднее, Σ_k – матрица ковариаций, |Σ_k| и $\Sigma _{k}^{{ - 1}}$ – определитель и обратная матрица матрицы ковариаций.

3. По пересчитанным вероятностям принадлежностей точек пересчитываются параметры распределений w_k, μ_k и Σ_k:

4. Пункты 2 и 3 алгоритма продолжаются до тех пор, пока алгоритм не сойдется, то есть параметры распределений не перестанут меняться.

Среди преимуществ метода GM можно отметить способность создания кластеров “эллиптической” формы, в отличие от методов KM и HW, в которых из-за минимизации внутрикластерного расстояния, при расчете которой дисперсии каждой из размерностей складываются в одно число (формула 2), кластеры получаются преимущественно “круглые”, т.е. кластеризация в них происходит по всем переменных независимо от их распределения.

В процессе EM-алгоритма происходит максимизация вероятности принадлежностей точек кластерам, а не минимизация V (формула 2), т.е. подбираются такие параметры K многомерных нормальных распределений, которые с наибольшей вероятностью описывают набор входных данных, что позволяет получать кластеры, имеющие вытянутую (“эллиптическую”) форму вдоль тех размерностей пространства набора данных, в которых выше дисперсия.

При режимном подходе к анализу атмосферной циркуляции метод кластеризации GM используется, например, в работах [Smyth et al., 1999; Hannachi, 1997], где с помощью него выделяются режимы низкочастотной изменчивости атмосферной циркуляции в тихоокеанском и атлантическом секторах Северного полушария.

4) Метод самоорганизующихся карт

Самоорганизующиеся карты Кохонена (Self-organizing maps, SOM) [Kohonen, 2012] представляют собой метод кластерного анализа, основанный на нейронных сетях. Данный метод позволяет получать кластеры из многомерного набора данных, которые удобно визуализировать на двумерной карте, кластеры в которой отсортированы друг относительно друга (отсюда “самоорганизующиеся”) на основе выбранной метрики расстояния (обычно Евклидовой).

Так как SOM входит в раздел методов машинного обучения на основе нейронных сетей, для него имеется своя сложившаяся терминология – кластеры называются “нейронами”. Перед запуском алгоритма кластеризации нейроны задаются в своем собственном двумерном пространстве, независимом от набора данных, т.е. на нейронной сети или карте. На нейронной сети определяются расстояния между нейронами или, иначе говоря, мера соседства нейронов, которая используется в процессе кластеризации для их взаимного “обучения”. Центроиды, т.е. координаты центров будущих кластеров, называются “векторами веса” или просто “весами” нейронов, итерации алгоритма называются “эпохами”.

Нейроны представляют из себя объекты, состоящие из двух векторов, заданных в разных пространствах – вектора веса, заданного в пространстве набора данных, и координат (как правило, двумерных), заданных в отдельном пространстве и определяющих положение нейрона на карте нейронной сети относительно других нейронов. В процессе работы алгоритма векторы веса нейронов как бы “натягиваются” на “облако” набора входных данных, при этом векторы веса нейронов, расположенных рядом друг с другом на карте нейронной сети, “притягиваются”, т.е. обновляют значения векторов веса относительно входных данных, к соседним областям в пространстве набора данных.

Алгоритм SOM состоит из следующих шагов:

1. Инициализация координат и весов нейронов. Координаты задаются на двумерной нейронной сети размера K = M × N, где M и N – размеры сети по горизонтали и вертикали. Нейроны могут взаимно располагаться в виде квадратной, шестиугольной или случайной сетки, от чего будет зависеть расстояние между соседствующими нейронами. Веса нейронов задаются (инициализируются) некоторым случайным образом в диапазоне значений набора входных данных, то есть в диапазоне координат, которые принимают точки в пространстве входных данных.

2. Из набора данных выбирается случайная точка. Среди кластеров (нейронов) ищется так называемый “нейрон-победитель”, вес которого на основе выбранной метрики, как правило, Евклидовой, будет ближе всего к данной точке:

где X_i – случайно выбранная точка из набора данных, W_c(t) – вес нейрона-победителя на эпохе t, W_k(t) – веса остальных нейронов.

3. Веса нейронов обновляются на основе расстояния от их весов до точки X_i и функции обучения, зависящей от эпохи и взаимного расстояния нейронов:

где W_k(t + 1) – обновленный вес k-го нейрона; (X_i – W_k(t)) – вектор (а не Евклидово расстояние), направленный от веса нейрона в сторону выбранной точки; h_ck(t) – функция обучения; ||r_c – r_k|| – расстояние между нейронами в пространстве нейронной сети, r – их двумерные координаты; α(t) и σ(t) – обучающие коэффициенты, монотонно убывающие с номером эпохи с тем, чтобы в начале работы алгоритма веса нейронов быстро обновлялись, а веса соседствующих нейронов притягивались к соседним областям пространства входных данных, но по мере работы алгоритма взаимное “обучение” затухало, и алгоритм сходился.

4. Пункты 2 и 3 алгоритма продолжаются до достижения некоторого критерия – либо пока не будет достигнут заранее предзаданный номер эпохи t_max (наиболее частый критерий), либо пока “ошибка карты” не перестанет существенно убывать: ΔV = V(t + 1) – V(t) ≤ ε; V можно посчитать по формуле 2, где в качестве C_k используются веса нейронов W_k(t).

Примеры использования нейронных сетей для выделения режимов крупномасштабной атмосферной циркуляции можно найти в работах [Polo et al., 2011; Loikith et al., 2017]. В работе [Polo et al., 2011] утверждается, что KM и SOM показывают похожие результаты, что также является одним из выводов настоящей работой (см. главу “результаты”).

Одно из основных преимуществ SOM состоит в автоматической сортировке кластеров друг относительно друга, то есть близкие по Евклидовому расстоянию кластеры будут располагаться рядом на нейронной карте. Например, на двумерной нейронной карте M × N = 5 × 5 (итого, K = 25 нейронов, то есть кластеров), кластеры с координатами {1,1}, {1,2} и {2,1} будут похожи друг на друга и иметь схожие характеристики, а кластеры с координатами {1,1} и {5,5} будут наиболее различаться между собой (по метрике Евклидова расстояния). Для методов KM и HW такая взаимная сортировка кластеров потребовала бы отдельной процедуры, то есть при кластеризации с таким же K = 25 кластеры с номерами 1 и 2 могли бы располагаться где угодно (в пространстве набора данных) относительно друг друга. При использовании SOM обычно выбирают большие K, т.к. для небольших K, например <10, взаимная сортировка, как правило, не нужна.

1) Метод локтя (Elbow method)

Данный метод является самым простым. Данные кластеризуются для каждого K в некотором интервале от 2 до K_max. Строится график функции V(K) (формула 2), для которой ищется “локоть” – точка перегиба, такое K, после которого V(K) перестает заметно убывать, т.е. последующее добавление новых кластеров не приводит к заметному улучшению качества кластеризации. Стоит отметить, что V(K) это, в любом случае, монотонно убывающая функция, которая стремится к 0 при K стремящемся к N.

2) Метод Силуэта

Для каждого потенциального K проводится кластеризация данных. Для каждой точки i из набора данных считается величина S(i), называемая коэффициент силуэта:

где N_k – число точек в кластере k, т.е. a(i) – среднее расстояние от точки до остальных точек из того же кластера, b(i) – среднее расстояние от точки из кластера k до точек из ближайшего для нее соседнего кластера k', S(i) – это величина, меняющаяся в пределах от –1 до 1, причем чем ближе она к 1, тем ближе точка к центроиду своего кластера и дальше от точек из других кластеров. Посчитав среднее или сумму S(i) по i для каждой из кластеризаций по всем потенциальным K, можем построить график от K. В отличие от метода локтя, при котором ищется точка перегиба функции V(k), в данном случае на графике ищется максимум <S(i)> от K.

3) Индекс классифицируемости

Индекс предложен в работе [Michelangeli et al., 1995] для выделения погодных режимов посредством кластеризации полей высоты геопотенциала на уровне 700 мбар. Идея состоит в том, что при оптимальном K различные решения при кластеризации методом KM, полученные путем разных инициализаций, более похожи друг на друга, чем при других K.

Допустим, что для данного K методом k-means по разным инициализациям посчитано M различных решений. Для каждой пары решений P, Q можно посчитать матрицу попарных корреляций центроидов кластеров:

где C – один из K центроидов данного решения (формула 2). Далее, для каждой строки i матрицы A_ij ищется максимум $A_{i}^{'}(P,Q)$ = $\mathop {\max }\limits_j {{A}_{{ij}}}(P,Q)$, т.е. $A_{i}^{'}(P,Q)$ – это лучшие попарные корреляции для каждого из центроидов кластеров в решении P, найденные среди центроидов кластеров в решении Q, т.е. корреляции центроидов похожих, аналогичных кластеров. Далее, среди пар аналогичных центроидов ищется “худшая” пара аналогичных центроидов, c(P, Q) = $\mathop {\min }\limits_i A_{i}^{'}(P,Q)$ = $\mathop {\min }\limits_i \left( {\mathop {\max }\limits_j {{A}_{{ij}}}(P,Q)} \right)$, которая и будет характеризовать схожесть двух различных решений (P и Q, в данном случае). Стоит отметить, что c(P, Q) ≠ c(Q, P). Посчитав среднее $c({{P}_{m}}(K),{{P}_{{m{\kern 1pt} '}}}(K))$ для всех имеющихся пар решений, получим индекс классифицируемости:

Для оценки значимости, полученный индекс классифицируемости c * (K) сравнивают с таковым, полученным по случайным сгенерированным выборкам. Выборки генерируются как марковские процессы, имеющие такую же автокорреляцию 1-го порядка и ковариацию, что и оригинальные данные. Полученные N_ген индексов с_ген * (K, N) сортируются по возрастанию. Долю с_ген * (K, N), которые оказались меньше с * (K), можно считать уровнем значимости индекса классифицируемости. Выбирается такое наименьшее K, для которого уровень значимости оказался выше уровня, заданного исследователем. При использовании данного индекса рекомендуется использовать как можно большие M и N_ген (50 и 100 в оригинальной статье).

4) Индекс воспроизводимости

Идея данного индекса, также описанного в [Michelangeli et al., 1995], заключается в том, что если для определенного K набор данных хорошо кластеризуется, то и его случайные половинчатые подвыборки будут сходиться к решениям, аналогичным таковым у оригинальной выборки.

Методом KM создается M различных решений, полученных с помощью разных инициализаций. Среди них ищется эталонное решение – такое m решение, для которого будет максимально среднее ❬c(P_m, P_m')❭_{m' ≠m} по всем остальным решениям m', т.е. решение, центроиды которого в среднем больше всего похожи на аналогичные центроиды остальных решений.

Далее, создаются R подвыборок, содержащих по 50% случайно выбранных точек оригинальной выборки данных. Для каждой из подвыборок, таким же образом, как для оригинальной выборки, проводится M кластеризаций и аналогичным образом выбирается эталонное решение. Методом попарных корреляций центроидов кластеров, аналогично индексу классифицируемости, для каждой пары эталонных решений оригинальной выборки P_m и одной из подвыборок P_mr считается $A_{i}^{'}({{P}_{m}},{{P}_{{mr}}})$ = $\mathop {\max }\limits_j {{A}_{{ij}}}({{P}_{m}},{{P}_{{mr}}})$ и осредняется по всем r = 1…R. Данная процедура проводится для каждого рассматриваемого K. Оптимальным можно считать такое K, для которого средняя по всем подвыборкам попарная корреляция центроидов эталонного решения оригинальной выборки и эталонных решений подвыборок – ${{\left\langle {A_{i}^{'}({{P}_{m}},{{P}_{{mr}}})} \right\rangle }_{r}}$ – будет выше некоторого заданного уровня значимости для каждого i = 1…K центроида эталонного решения оригинальной выборки.

Вопрос количества погодных режимов, остается одним из ключевых при режимном подходе в исследованиях крупномасштабной атмосферной циркуляции [Christiansen, 2007; Philipp et al., 2007]. Оптимальное K, если оно существует, зависит как от выбранного метода кластеризации, так и от свойств набора данных – исследуемой характеристики, определяющей циркуляцию (поле давления на уровне моря, высоты геопотенциальной поверхности или другой переменной), выбранной области, продолжительности временного ряда, процедуры предобработки данных и т.д. [Falkena et al., 2020].

1) Данные и предобработка

Для идентификации режимов использовались поля высоты геопотенциала на уровне 500 гПа (z500) реанализа ERA5 [Hersbach et al., 2020] c 1940 по 2022 гг. Поля были осреднены до суточных и интерполированы на сетку с шагом в 1 градус для удобства. Для выделения режимов в Евро-Атлантике был выбран сектор 80 з.д.–40 в.д. по долготе и 30 с.ш.–80 с.ш. по широте, такая область чаще всего используется для выделения погодных режимов в ЕАТ [Fabiano et al., 2020].

Перед кластеризацией были получены поля аномалий путем вычитания из суточных полей сезонного хода, который считался осреднением суточных полей за одни и те же календарные даты разных лет за весь период (1940–2022 гг.), сглаженные скользящим средним в 5 суток. Для удаления высокочастотных синоптических волн к полям аномалий применялся фильтр Баттерворта низких частот [Selesnick et al., 1998] с периодом фильтрации 10 суток. Из фильтрованных суточных полей аномалий высоты геопотенциала выбирались зимние (декабрь, январь, февраль) и летние (июнь, июль, август) месяцы для кластеризации зимних и летних режимов соответственно. За исследуемый 83-летний период 1940–2022 гг. для зимних месяцев имеется 7491 суточных полей, для летних – 7636.

Для сглаживания пространственной дисперсии полей и уменьшения размерности данных фильтрованные суточные поля аномалий высоты геопотенциала (далее, поля аномалий) раскладывались на эмпирические ортогональные функции (ЭОФ). Перед разложением на ЭОФ поля аномалий были взвешены на корень из площади для уравнивания вклада в дисперсию узлов сетки полей на разных широтах. Для зимних полей было отобрано 13 первых ЭОФ, объясняющих 94.6% дисперсии, для летних полей 20 первых ЭОФ, объясняющих 94.8% дисперсии, такие числа ЭОФ выбраны, чтобы они объясняли примерно по 95% дисперсии в обоих случаях. Непосредственно кластеризации подвергались временные ряды главных компонент ЭОФ.

Отметим, что в разных работах с выделением режимов атмосферной циркуляции в ЕАТ перед кластеризацией используют разное число ЭОФ, например, 4 в работах [Dawson et al., 2012; Fabiano et al., 2020] или 14 в [Cassou, 2008; Charlton-Perez et al., 2018]. Как правило, авторы отмечают, что результаты кластеризации не сильно меняются при увеличении числа ЭОФ. В работе [Falkena et al., 2020] исследуются различия в погодных режимах, полученных методом KM, в зависимости от способа предобработки данных – использования (или неиспользования) временных фильтров и разложения на ЭОФ. Утверждается, что режимы, полученные при разложении данных на 10, 15 и 20 первых ЭОФ и без разложения (оригинальные данные), практически не отличаются и имеют схожие повторяемости и вероятности переходов.

2) Выбор количества режимов

Была проведена оценка оптимального количества режимов на интервале K от 2 до 10 для кластеризации методом KM. Как было отмечено выше, простые методы проверки оптимального числа кластеров (метод локтя и Силуэт) не дают четкого ответа. Графики V(K) и $\left\langle {{{S}_{i}}\left( t \right)} \right\rangle $ от K представляют собой убывающие функции, однако для зимних режимов наблюдается слабый перегиб V(K) и локальный максимум $\left\langle {{{S}_{i}}\left( t \right)} \right\rangle $ (при данном методе ищется глобальный максимум) при k = 4. Для летних режимов аналогичной картины не наблюдается, кривые V(K) и $\left\langle {{{S}_{i}}\left( t \right)} \right\rangle $ от K плавно убывают.

Для оценки оптимального числа режимов использовались в том числе индексы классифицируемости и воспроизводимости. Для подсчета индекса классифицируемости было сгенерировано 400 случайных выборок с одинаковыми дисперсиями и автокорреляцией 1-го порядка, как у оригинальных данных, и 400 случайных подвыборок для подсчета индекса воспроизводимости.

Для зимних режимов индексы показывают локальные пики для k = 4 и k = 6 (см. рис. 1), но ни для одного K индекс классифицируемости оригинальной выборки не оказался выше, чем у 95% сгенерированных выборок, т.е. никакое K не будет являться статистически значимым. Для летних режимов высокое значение индексов наблюдается при k = 3 (см. рис. 2).

В целом, использованные методы не дали однозначного ответа об оптимальном количестве режимов. Так как в большинстве исследований для выделения зимних режимов в Евро-Атлантике используется K = 4 [Cassou, 2008; Michelangeli et al., 1995; Fabiano et al., 2020; Dawson et al., 2012; Charlton-Perez et al., 2018], было решено выбрать это число для возможности сопоставления с предыдущими результатами. Для летних режимов кластеризация так же проводилась с K = 4 для возможности сравнения полученных летних режимов с зимними. Стоит отметить, что существуют работы, в которых в ЕАТ зимой выделяется 6 режимов [Falkena et al., 2020], в таком случае режимы “NAO+”, “NAO–” остаются визуально неизменными (поля центроидов этих режимов), а у режимов “SB” и “AR” появляются “противоположные” – “SB-” и “AR-”.

Обработка данных и их последующая кластеризация выполнялись в программной среде MATLAB с использованием как встроенного функционала, так и дополнительных программ, разработанных в том числе авторами статьи. Предобработка данных (удаление сезонного хода, фильтрация по времени, разложение на ЭОФ) выполнялась как с использованием встроенных функций, так и функций из пакета “Climate Data Toolbox for MATLAB” [Greene et al., 2019] (“filt1”, “eof” и другие). Для кластеризации KM с использованием оптимизаций “kmeans++” и “simulated annealing”, а также для подсчета индексов классифицируемости и воспроизводимости использовались написанные авторами программы, для остальных методов кластеризации – GM, HW и SOM – использовались встроенные в MATLAB функции (“fitgmdist”, “cluster”, “linkage”, “nctool” и прочие).

1) Поля зимних и летних погодных режимов

Средние характерные поля аномалий z500 режимов, полученные путем усреднения всех cуточных полей аномалий z500, попавших в соответствующие режимы, представлены на рис. 3. Режимы, полученные методом SOM, не приведены на рисунке. Результаты кластеризации методом SOM оказались практически идентичны результатам, полученным по методу KM – число суточных полей, попавших в другие, чем в методе KM, режимы, составило 12 и 10 полей из 7491 и 7636 для зимних и летних режимов соответственно, то есть менее 0.2% (также см. табл. 1), поэтому далее в данной главе результаты метода SOM не приведены, так как принимаются идентичными таковым по методу KM.

Традиционные названия NAO+, NAO–, SB, AR даны для режимов, полученных методом KM, так как визуально они наиболее совпадают с традиционными зимними ЕАТ режимами, полученными этим методом в других исследованиях [Fabiano et al., 2020; Dawson et al., 2012; Charlton-Perez et al., 2018]. Не все режимы, полученные по методам HW и GM, воспроизводят вышеназванные режимы, однако они сопоставлены с режимами по методу KM таким образом, чтобы средняя пространственная корреляция режимов, названных аналогично, была максимальной. На рис. 3 и далее по тексту для летних режимов используются аналогичные зимним режимам (NAO+, NAO–, SB, AR) названия с припиской “s” (summer).

Из рис. 3 видно, что летние поля режимов являются менее выраженными (поля на рисунках более бледные даже с учетом использования уменьшенной шкалы) по сравнению с зимними, что является следствием более низкой изменчивости полей z500 в летние месяцы по сравнению с зимними.

В табл. 1 приведено сравнение внутрикластерных расстояний в зависимости от выбранного метода кластеризации. Видно, что KM и SOM произвели кластеры с близкими внутрикластерными расстояниями (разница в пользу k-means в 5-ом знаке после запятой, не приведено в таблице), при этом кластеризации HW и GM оказываются хуже в терминах суммарных внутрикластерных расстояний и доли объясненной дисперсии, то есть в среднем кластеры HW и GM менее плотные, несмотря на выраженность некоторых из них (например, NAO– в GM). Стоит отметить, что доля объясненной дисперсии (которая считается как 1 – V_k/V₁) зависит от метода предобработки данных и выбранного числа ЭОФ разложения оригинальных полей, однако тенденция превосходства k-means при этом сохраняется.

В табл. 2 представлено сравнение пространственных корреляций средних полей режимов, полученных разными методами, а также сравнение пространственных корреляций зимних и летних полей для одного и того же метода. По приведенным значениям в таблице видно, что результаты кластеризации GM и HW плохо коррелируют между собой (<0.5 в среднем для обоих сезонов), но хорошо коррелируют с KM (>0.7 в среднем для обоих сезонов). Сравнение пространственных корреляций зимних и летних полей аналогичных режимов показывает, что наиболее похожие летние режимы по сравнению с зимними получаются при методе KM, хуже всего при методе HW, что так же свидетельствует о том, что KM является более предпочтительным методом при кластеризации полей циркуляции атмосферы.

2) Относительная и сезонная повторяемость режимов

Для полученных разными методами режимов анализировались и сравнивались относительная и сезонная повторяемости, средняя продолжительность и матрицы переходов. Относительная повторяемость режимов представлена в табл. 3, где также приведена повторяемость классических зимних погодных режимов, полученных другими авторами кластеризацией KM суточных полей z500 за разные временные периоды и с немного отличающимися способами предобработки данных.

По табл. 3 видно, что полученная авторами повторяемость режимов NAO+/NAO– в пределах 1.5% согласуется с таковой в других работах [Cassou, 2008; Dawson et al., 2012; Charlton-Perez et al., 2018]. Полученные авторами повторяемости режимов SB и AR сильнее отличаются от их повторяемостей в других работах, что, скорее всего, является следствием большей чувствительности этих режимов к выбору временного интервала (1940–2022 гг. в настоящей работе) и способу предобработки данных (13 ЭОФ, 10-дневный фильтр). Прослеживается общая тенденция – самый частый режим – NAO+, средние по частоте SB и AR, NAO– – самый редкий. При кластеризации z500 другими методами данная тенденция не прослеживается – в HW режим NAO– оказался наиболее продолжительным с частотой в 33.38%, в GM NAO– – наименее продолжительный, но со слишком низкой относительной повторяемостью 12.25%.

Несмотря на хорошую пространственную корреляцию средних полей летних режимов по сравнению с зимними при методе KM (средняя корреляция 0.84, см. табл. 2), их относительная повторяемость отличается – самым частым оказался режим sNAO– с повторяемостью 28.06% против 21.39% зимой, а повторяемость sNAO+ уменьшилась с 30.20 до 25.51%. В целом, относительная населенность летних режимов получилась более равномерной по сравнению с зимними, что является одним из свойств метода KM при плохой кластеризуемости данных (доля объясненной кластеризацией дисперсии 23.62% летом против 32.64% зимой, см. табл. 1). Для методов HW и GM отмечается излишняя населенность некоторых режимов в летние месяцы – sNAO– для HW и sNAO+ для GM.

Помимо относительной повторяемости режимов, анализировались временные ряды сезонной повторяемости (количество суток за сезон, когда наблюдался данный режим) и их линейные тренды. Для сезонной повторяемости зимних режимов, полученных разными методами, получились разные по величине и значимости тренды, но наблюдается общая картина – положительный тренд повторяемости NAO+, который оказался значимым (на уровне 95%) для режимов, полученных по методам KM и GM, и отрицательный тренд повторяемости NAO–, который оказался значимым для режимов, полученных по методам KM и HW.

Временные ряды и линейные тренды режимов NAO+ и NAO–, полученных методом KM, представлены на рис. 4. Режимы NAO+ и NAO– ассоциированы с волнами тепла и холода соответственно над территориями Европы и Западной России [Cattiaux et al., 2010; Бардин и др., 2019], поэтому их многолетние тренды сезонной повторяемости могут частично объяснять общий тренд потепления в данных регионах в последние десятилетия.

Зимние режимы SB и AR, полученные по различным методам, имеют незначимые разнонаправленные тренды сезонной повторяемости, кроме значимого отрицательного тренда SB при методе GM, который, однако, слабо похож на свои аналоги в методах HW и KM (см. рис. 3 и табл. 2), т.е., вообще говоря, не представляет собой скандинавский блокинг. Для летних месяцев значимые тренды сезонной повторяемости отсутствуют практически для всех режимов, полученных разными методами, кроме положительного тренда для режима sNAO–, полученного по методу GM.

3) Характерная продолжительность режимов

В табл. 4 представлены средние продолжительности режимов. Самым продолжительным оказался режим NAO–, причем как для зимних месяцев (>10 дней по методам KM и HW), так и для летних (sNAO–, >9 дней по методам KM и HW). Такая повышенная продолжительность отмечается на фоне самой низкой (по методу KM) среди четырех режимов относительной повторяемости в зимние месяцы. Наименее продолжительными оказались режимы SB в зимние месяцы и sAR в летние.

Сравнение средних и медианных значений продолжительности режимов, полученных разными методами, показывает, что наиболее продолжительные режимы получаются при кластеризации методом HW, а наименее продолжительные при кластеризации GM. Для всех режимов отмечается высокая величина стандартного отклонения их продолжительности, которая близка к средней продолжительности и превышает медианную, что говорит о том, что среди событий режимов часто встречаются как и очень короткие, продолжительностью не более трех суток, которые можно условно назвать переходными, так и довольно продолжительные, которые могут занимать до двух месяцев (средняя максимальная продолжительность по всем режимам и методам для зимних месяцев составляет 45 сут, для летних 40 сут). На рис. 5 представлены гистограммы распределения продолжительностей для режимов, полученных по методу KM. Из рисунка 5 и табл. 4 видно, что летние режимы оказываются менее продолжительными, чем зимние, также по рисунку заметно отличие в распределении продолжительностей режима NAO– и sNAO– относительно остальных.

Стоит упомянуть, что характерная продолжительность режимов сильно зависит от способа обработки данных, а именно от применения или неприменения фильтрации по времени и выбора периода отсечки. В данной работе применялся фильтр Баттерворта низких частот с периодом отсечки в 10 сут, из-за которого общая дисперсия суточных полей аномалий z500 снижается примерно в 1.5 раза, а характерная продолжительность режимов увеличивается.

4) Переходы между режимами

Были рассчитаны матрицы переходов между режимами, а также статистически значимые переходы (см. главу 4) с целью выяснить, какие переходы более или менее вероятны, что может иметь потенциальное применение в задачах предсказуемости. Вероятностные матрицы переходов режимов, в том числе статистически значимые из них, полученные для каждого из методов, представлены в табл. 5. Вероятностные матрицы получены путем деления числа переходов из данного режима в другие на суммарное число переходов из данного режима (т.е. сумма вдоль каждой строки в табл. 5 равна 1 без учета особенностей округления).

Статистические значимые переходы получены по методике, предложенной в [Vautard et al., 1990], т.е. это такие переходы, которые происходят в 95% случаев чаще или реже, чем в матрицах, случайно сгенерированных путем перемешивания режимных событий с сохранением общего числа событий каждого режима. Так как переходы рассматривались между режимными событиями, а не между отдельными суточными полями, то переходы режимов в самих себя (например, NAO+ → NAO+) при таком подходе отсутствуют.

Из табл. 5 видно, что вероятность перехода из режима NAO+ в “противоположный” режим NAO- довольно низкая – 0.11 для режимов по методу KM и 0.12 по GM, однако 0.36 по HW. При этом вероятность перехода из NAO– в NAO+ в 3–4 раза выше, чем из NAO+ в NAO–, что может объясняться как свойствами атмосферной циркуляции в Евро-Атлантике, так и разницей в повторяемости этих режимов, так как более населенные режимы “перетягивают” переходы в свою сторону, однако NAO+ является более населенным только по методам KM и GM, из-за чего, вероятно, такой разницы в вероятности переходов не наблюдается для этих режимов, полученных методом HW (0.36 для перехода из NAO+ в NAO– и 0.35 для перехода из NAO+ в NAO–), при этом режим NAO+ в HW заметно отличается от NAO+ по методам KM и GM (см. табл. 2).

Аналогичный, хоть и менее выраженный результат получился для летних режимов sNAO+ и sNAO– – вероятность перехода из sNAO– в sNAO+ 1.5–2 раза выше (для каждого из методов), чем из sNAO+ в sNAO–, что, при схожести летних и зимних режимов NAO разных фаз, может свидетельствовать о том, что разница в вероятности переходов между этими режимами объясняется не только разницей в их повторяемости, так как в летние месяцы повторяемость sNAO– выше, чем sNAO+ для всех методов, кроме GM, о чем так же свидетельствует тот факт, что для всех методов переход из sNAO– в sNAO+ оказался статистически значимо частым, а переход из sNAO+ в sNAO– статистически значимо редким.

Интересным результатом оказалась нулевая вероятность перехода между режимами NAO– и AR при кластеризации методом GM, т.е. за исследуемый период между полученными режимами такие переходы не наблюдались ни разу, однако стоит отметить, что режим AR для метода GM отличен от классического режима AR, полученного с помощью KM (см. рис. 3 и табл. 2). При других методах кластеризации между этими режимами наблюдаются либо незначимые, либо статистически значимо частые (NAO– → AR в HW) вероятности переходов. Также для метода GM, в отличие от остальных, наблюдаются “цепочки” значимых переходов, проходящие через каждый из режимов: NAO+ → AR → SB → NAO– → NAO+ → →  … и, аналогично, sNAO+ → sAR → sSB → → sNAO– → sNAO+ → …, т.е. системы режимов, полученных методов GM для зимних и летних месяцев, со статистически значимо высокой вероятностью поочередно проходят через все свои состояния, несмотря на то, что летние режимы при методе GM относительно плохо воспроизвелись по сравнению с зимними (средняя корреляция полей аналогичных режимов 0.6). Такой результат может являться как простым совпадением, так и свойством метода GM.

Для других методов не наблюдается замкнутых цепочек статистически значимо частых переходов, проходящих через все режимы, так как не для каждого из режимов, полученных по методам KM и HW существуют значимо частые переходы, например, их нет для режимов NAO– и sAR по методу KM и режимов AR и sNAO+ по методу HW.