Журнал вычислительной математики и математической физики, 2021, T. 61, № 7, стр. 1113-1124

1. ВВЕДЕНИЕ

В настоящее время обработка и анализ результатов длительных измерений, проводимых через заданные промежутки времени, вызывают интерес в различных сферах, от экономики до метеорологии. Для этого используется достаточно обширный арсенал методов как математически обоснованных (например, метод Фурье для стационарных рядов, вейвлет-преобразования для нестационарных), так и эвристических [1]. В этой ситуации интерес представляет разработка методов анализа временных рядов, основанных на достаточно общих математических моделях, дающих оптимальные решения на широком классе задач.

В данной работе для исследования цикличности временных рядов в первую очередь используется морфологический подход, развиваемый в школе Ю.П. Пытьева [2]–[4]. Предлагается способ морфологической фильтрации, позволяющий выделить в сигнале составляющую, которая представляет собой последовательность фрагментов, схожих по форме, но меняющихся по длительности и амплитуде в некоторых заданных пределах. Такая фильтрация может быть полезна для декомпозиции исходного временного ряда на “почти циклическую” составляющую, связанную с сезонностью, и остаток, который может быть подвергнут дальнейшему исследованию с целью выделения дополнительных циклических составляющих.

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

2. МОРФОЛОГИЧЕСКИЙ ПОДХОД К АНАЛИЗУ СИГНАЛОВ

На практике, в частности, при дистанционном зондировании атмосферы, при изучении временных рядов метеорологических параметров и др., анализируемый сигнал имеет вид суммы почти повторяющихся фрагментов, схожих по форме с сигналами иной формы – со стационарным или нестационарным сигналами, с широкополосной шумовой составляющей и др. Задача состоит в разделении суммы этих сигналов на составляющие и анализе параметров каждой из них.

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

Более гибкую модель сигнала заданной формы дает морфологический подход [2]–[4]. Морфологические методы анализа сигналов, развиваемые в школе Ю.П. Пытьева, предназначены для анализа структуры сигнала, сохраняющейся при его преобразовании, принадлежащем заданному классу. В частности, если каждый из двух сигналов представляет собой набор “пиков” и “провалов” (т.е. локальных максимумов и минимумов), то считать их сигналами одной формы разумно, если положения “пиков” и “провалов” совпадают, а упорядоченность локальных экстремумов по амплитуде сохраняется (самый высокий “пик” у обoих сигналов расположен в одной и той же точке и т.д.). Это свойство сигналов сохранится при любом монотонном преобразовании их амплитуд.

Похожий подход к анализу сигналов и изображений предложен в работах А.Н. Каркищенко и A.B. Гончарова [6] и [7], используемое в них описание изображений и сигналов получило название “знакового описания”. Изучению свойств такого описания посвящена монография [8], исследования продолжены в работе [9]. Методы, в которых используется свойство неотрицательности производных сигнала, развиваются также в работе [10].

В методах морфологического анализа [3] форма сигнала $f(t)$, $t \in T$, понимается как множество ${{V}_{f}}$ сигналов, полученных из $f(t)$ всевозможными монотонными преобразованиями его амплитуды; очевидно, при таком понимании формы указанное выше характерное свойство сигнала сохраняется. Математическую модель сигнала и множество преобразований, сохраняющих форму, обычно выбирают так, чтобы для любого предъявленного для анализа сигнала $g(t)$, $t \in T$, была разрешима задача его наилучшего приближения элементами множества ${{V}_{f}}$, решение которой называют проекцией $g$ на ${{V}_{f}}$. Операцию проецирования, в ряде случаев определенную конструктивно, называют формой сигнала (изображения) $g$ [3].

В настоящей работе методы морфологического анализа, описанные в [2]–[4], адаптируются для анализа временного ряда концентрации СО₂, полученного на основе данных круглогодичных измерений на эколого-климатической станции “AsiaFlux” во Вьетнаме с 2011 по 2018 г., предоставленных ИПЭЭ им. А.Н. Северцова РАН [11]. Концентрация CO₂ в атмосфере регистрируется ежесекундно с последующим усреднением до получасовых значений и выражена в миллионных долях (ppm). Таким образом, рассматриваемый временной ряд концентрации CO₂ представляет собой упорядоченную последовательность усредненных получасовых значений концентрации CO₂ на определенной высоте над поверхностью Земли. Во временном ряду имеется составляющая, отражающая циклические процессы с периодом в одни сутки, амплитуда этой составляющей существенно больше других составляющих временного ряда. При этом от суток к суткам имеется существенная вариация поведения концентрации СО₂, что затрудняет применение не только фурье-, но и вейвлет-анализа. В данной работе предлагается математическая модель формы сигнала, отражающая суточную динамику концентрации CO₂, в виде функции, направление выпуклости которой меняется на противоположное в точках перегиба. Эти точки являются параметрами формы фрагмента сигнала, моделирующего суточный ход концентрации CO₂. Считается, что два фрагмента сигнала имеют одинаковую форму, если интервалы их выпуклостей вверх и выпуклостей вниз совпадают.

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

3. МАТЕМАТИЧЕСКАЯ МОДЕЛЬ ВРЕМЕННОГО РЯДА

Временной ряд будем рассматривать как конечный набор значений ${{\xi }_{i}}$, $i = 1,\; \ldots ,\;N$. Здесь каждое значение ${{\xi }_{i}}$ есть результат наблюдения некоторого параметра в i-й момент времени. Зависимость значения параметра от времени является суммой почти циклической составляющей ${{f}_{i}}$ и широкополосной случайной помехи ${{\nu }_{i}}$, $i = 1,\; \ldots ,\;N$.

Зададим математическую модель почти циклической составляющей ${{f}_{i}}$, $i = 1,\; \ldots ,\;N$, считая, что последовательность моментов {1, 2, …, N} времени измерений можно разбить на участки {${{i}_{1}},\; \ldots ,\;{{i}_{2}} - 1$},{${{i}_{2}},\; \ldots ,\;{{i}_{3}} - 1$}, …, {${{i}_{{2m - 1}}},\; \ldots ,\;.{{i}_{{2m}}} - 1$}, длительность которых может изменяться от участка к участку в заданных пределах, и на каждом нечетном участке почти периодическая составляющая имеет выпуклость вверх, а на четном – вниз. В точках перегиба при переходе от выпуклости вниз к выпуклости вверх функция не убывает, а при смене выпуклости вверх на выпуклость вниз – не возрастает. Множество ${{V}_{f}}({{i}_{1}},\; \ldots ,\;{{i}_{{2m}}})$ всех таких функций назовем формой почти периодического сигнала, а границы участков $1 \leqslant {{i}_{1}} \leqslant \; \ldots \; \leqslant {{i}_{{2m}}} \leqslant N + 1$ определяют параметры его формы. Итак, форма ${{V}_{f}}({{i}_{1}},\; \ldots ,\;{{i}_{{2m}}})$ почти периодического сигнала есть набор значений ${{f}_{i}}$, $i = {{i}_{1}},\; \ldots ,\;{{i}_{{2m}}} - 1$, для которых выполнены неравенства

Если рассматривать значения ${{\xi }_{i}}$, $i = {{i}_{1}},\; \ldots ,\;{{i}_{{2m}}} - 1$, временного ряда как координаты вектора конечномерного евклидова пространства, то множество ${{V}_{f}}({{i}_{1}},\; \ldots ,\;{{i}_{{2m}}})$ является выпуклым замкнутым множеством, и для любого вектора ${\mathbf{\xi }}$ с координатами ${{\xi }_{i}}$, $i = {{i}_{1}},\; \ldots ,\;{{i}_{{2m}}} - 1$, однозначно разрешима задача наилучшего приближения вектора ${\mathbf{\xi }}$ векторами из ${{V}_{f}}({{i}_{1}},\; \ldots ,\;{{i}_{{2m}}})$ при фиксированных параметрах формы. Ее решение $P({{i}_{1}},\; \ldots ,\;{{i}_{{2m}}}){\mathbf{\xi }}$ называется проекцией ${\mathbf{\xi }}$ на ${{V}_{f}}({{i}_{1}},\; \ldots ,\;{{i}_{{2m}}})$, а оператор проецирования $P({{i}_{1}},\; \ldots ,\;{{i}_{{2m}}})$ взаимно однозначно определяется множеством ${{V}_{f}}({{i}_{1}},\; \ldots ,\;{{i}_{{2m}}})$, и поэтому тоже называется формой сигнала ${\mathbf{f}}$ [3].

Сформулируем лемму, на основе которой удается найти решение задачи (2).

Лемма. Пусть в евклидовых пространствах ${{R}^{m}}$, ${{R}^{k}}$, ${{R}^{{m + k}}} = {{R}^{m}} \otimes {{R}^{k}}$, где ${\mathbf{g}} = ({\mathbf{t}},{\mathbf{s}})$, ${\mathbf{g}} \in {{R}^{{m + k}}}$, ${\mathbf{t}} \in {{R}^{m}}$, ${\mathbf{s}} \in {{R}^{k}}$ и $({{g}_{1}},\; \ldots ,\;{{g}_{m}},{{g}_{{m + 1}}},\; \ldots ,\;{{g}_{{m + k}}}) = ({{t}_{1}},\; \ldots ,\;{{t}_{m}},{{s}_{1}},\; \ldots ,\;{{s}_{k}})$, выпуклые замкнутые конусы ${{V}_{{m + k}}}$, ${{V}_{m}}$ и ${{V}_{k}}$ заданы системами линейных неравенств: ${{V}_{m}} = \{ {\mathbf{t}} \in {{R}^{m}}:A{\mathbf{t}} \geqslant {\mathbf{c}}\} $, ${{V}_{k}} = \{ {\mathbf{s}} \in {{R}^{k}}:B{\mathbf{s}} \geqslant {\mathbf{d}}\} $,

для некоторых матриц $A,\;B,\;C$ и векторов ${\mathbf{c,}}\;{\mathbf{d,}}\;{\mathbf{g}}$; ${{P}_{{m + k}}}$, ${{P}_{m}}$ и ${{P}_{k}}$ – проекторы на конусы ${{V}_{{m + k}}}$, ${{V}_{m}}$ и ${{V}_{k}}$ соответственно. Тогда для некоторого вектора ${\mathbf{\tilde {g}}} = ({\mathbf{\tilde {t}}},{\mathbf{\tilde {s}}}) \in {{R}^{{m + k}}}$ равенство ${{P}_{{m + k}}}({\mathbf{\tilde {t}}},{\mathbf{\tilde {s}}}) = ({{P}_{m}}{\mathbf{\tilde {t}}},{{P}_{k}}{\mathbf{\tilde {s}}})$ выполняется тогда и только тогда, когда выполнено неравенство $C({{P}_{m}}{\mathbf{\tilde {t}}},{{P}_{k}}{\mathbf{\tilde {s}}}) \geqslant {\mathbf{g}}$.

Эта лемма обобщает утверждения, приведенные в работах [3] и [4], и доказывается аналогично.

Лемма позволяет последовательно строить проекции фрагментов временного ряда на конус, определяемый неравенствами (1). На первом шаге строится проекция фрагмента ряда, состоящего из трех первых последовательных значений, на конус размерности 3. Далее, проекции фрагментов из 4, 5, … значений элементов ряда на соответствующие конусы размерности 4, 5 и т.д. строятся путем последовательного присоединения к проецируемому фрагменту по одному значению.

Действительно, пусть для определенности ${{i}_{1}} = 1$, и построена проекция ${{P}_{j}}({{f}_{1}},\; \ldots ,\;{{f}_{j}})$ на конус ${{V}_{j}}$, определяемый неравенствами (1), в которые входят только координаты $({{f}_{1}},\; \ldots ,\;{{f}_{j}})$. Тогда проекция ${{P}_{{j + 1}}}({{f}_{1}},\; \ldots ,\;{{f}_{{j + 1}}})$ на конус ${{V}_{{j + 1}}}$ строится последовательным построением пары проекций ${{P}_{{j - l + 1}}}({{f}_{1}},\; \ldots ,\;{{f}_{{j - l + 1}}})$, и ${{P}_{l}}({{f}_{{j - l + 2}}},\; \ldots ,\;{{f}_{{j + 1}}})$ вектора $({{f}_{1}},\; \ldots ,\;{{f}_{{j + 1}}})$ на конусы ${{V}_{{j - l - 1}}}$ и ${{V}_{{l + 2}}}$ (l = 1, 2, …) соответственно. Здесь первый конус ${{V}_{{j - l - 1}}}$ определяется неравенствами (1), в которые входят только координаты $({{f}_{1}},\; \ldots ,\;{{f}_{{j - l + 1}}}),$ а второй ${{V}_{{l + 2}}}$ неравенствами (1), в которые входят только координаты $({{f}_{{j - l}}},\; \ldots ,\;{{f}_{{j + 1}}})$, l = 1, 2, …, до тех пор, пока не будут выполнены неравенства (1), в которые входит хотя бы одна координата как из первого набора, так и из второго. Заметим, что как в первый, так и во второй наборы координат могут входить одна, две и более координат, в первых двух случаях могут отсутствовать неравенства, связывающие только координаты из соответствующего набора; в этой ситуации ограничения на эти координаты отсутствуют, а конус превращается в линейное подпространство.

Рассмотрим подробнее, как строится проекция ${{P}_{{j + 1}}}({{f}_{1}},\; \ldots ,\;{{f}_{{j + 1}}})$ при известной ${{P}_{j}}({{f}_{1}},\; \ldots ,\;{{f}_{j}})$. Пусть имеется два набора координат: $({{f}_{1}},\; \ldots ,\;{{f}_{j}})$ и ${{f}_{{j + 1}}}$. Если выполнены все неравенства из (1), в которые входят ${{f}_{{j + 1}}}$ и хотя бы одна координата из $({{f}_{{j - 1}}},{{f}_{j}})$, то искомая проекция построена и равна ${{P}_{{j + 1}}}({{f}_{1}},\; \ldots ,\;{{f}_{{j + 1}}}) = ({{P}_{j}}({{f}_{1}},\; \ldots ,\;{{f}_{j}}),{{f}_{{j + 1}}})$, так как выполнены условия леммы при $B = 0$, ${\mathbf{d}} = 0$.

Если эти неравенства не выполняются, следует построить две проекции: ${{P}_{{j - 1}}}({{f}_{1}},\; \ldots ,\;{{f}_{{j - 1}}})$ (она уже вычислена на предыдущих шагах) вектора $({{f}_{1}},\; \ldots ,\;{{f}_{{j - 1}}})$ на конус ${{V}_{{j - l}}}$ и ${{P}_{2}}({{f}_{j}},{{f}_{{j + 1}}})$ вектора $({{f}_{j}},{{f}_{{j + 1}}})$ на множество, определяемое неравенствами из (1), в которые входят только координаты $({{f}_{j}},{{f}_{{j + 1}}})$; если таких неравенств нет, то ${{P}_{2}}({{f}_{j}},{{f}_{{j + 1}}}) = ({{f}_{j}},{{f}_{{j + 1}}})$. Далее для координат построенных проекций проверяются выполнения неравенств (1), в которые входят как координаты из набора $({{f}_{j}},{{f}_{{j + 1}}})$, так и из набора $({{f}_{{j - 2}}},{{f}_{{j - 1}}})$. Если эти неравенства выполняются, то искомая проекция построена, иначе следует продолжить процесс.

Поясним, как построить проекции на конусы на начальных шагах. Пусть, например, как и прежде, ${{i}_{1}} = 1$, а ${{i}_{2}} > 4$.

Шаг 1. Выберем первые три координаты ${{f}_{1}}$, ${{f}_{2}}$, ${{f}_{3}}$, и поcтроим проекцию ${{P}_{3}}({{f}_{1}},{{f}_{2}},{{f}_{3}}) = ({{q}_{1}},{{q}_{2}},{{q}_{3}})$ трехмерного вектора (${{f}_{1}}$, ${{f}_{2}}$, ${{f}_{3}}$) на трехмерный конус

решая задачу выпуклого программирования Эта проекция есть либо сам вектор (${{f}_{1}}$, ${{f}_{2}}$, ${{f}_{3}}$), если $2{{f}_{2}} \geqslant {{f}_{1}} + {{f}_{3}}$, либо наилучшая среднеквадратичная аппроксимация вектора (${{f}_{1}}$, ${{f}_{2}}$, ${{f}_{3}}$) вектором с координатами $({{q}_{1}},({{q}_{1}}, + {{q}_{3}}){\text{/}}2,{{q}_{3}})$, т.е. линейно зависящими от времени. Заметим, что эта ситуация соответствует утверждению леммы при $B = 0$, $C = 0$, ${\mathbf{d}} = 0$, ${\mathbf{g}} = 0$.

Шаг 2. Добавим к набору координат (${{f}_{1}}$, ${{f}_{2}}$, ${{f}_{3}}$) вектора ${\mathbf{f}}$ четвертую координату ${{f}_{4}}$. Если  выполнены неравенства $2{{q}_{3}} \geqslant {{q}_{2}} + {{f}_{4}}$ для координат проекции ${{P}_{3}}({{f}_{1}},{{f}_{2}},{{f}_{3}}) = ({{q}_{1}},{{q}_{2}},{{q}_{3}})$, вычисленной  на  первом  шаге,  и  ${{f}_{4}}$,  то  согласно  лемме  проекция  ${{P}_{{3 + 1}}}({{f}_{1}},{{f}_{2}},{{f}_{3}},{{f}_{4}})$ на конус ${{V}_{{3 + 1}}} = \{ ({{t}_{1}},{{t}_{2}},{{t}_{3}},{{s}_{1}}):2{{t}_{2}} \geqslant {{t}_{1}} + {{t}_{3}},\;2{{t}_{3}} \geqslant {{t}_{2}} + {{s}_{1}}\} $ в четырехмерном пространстве равна ${{P}_{{3 + 1}}}({{f}_{1}},{{f}_{2}},{{f}_{3}},{{f}_{4}}) = ({{P}_{3}}({{f}_{1}},{{f}_{2}},{{f}_{3}}){\mathbf{,}}{{f}_{4}})$ $ = ({{P}_{3}}({{f}_{1}},{{f}_{2}},{{f}_{3}}),{{f}_{4}})$, что соответствует условиям леммы при ${\mathbf{d}} = 0$. Если это неравенство не выполнено, то дальнейшие действия зависят от результата первого шага.

1. Если ${{P}_{3}}({{f}_{1}},{{f}_{2}},{{f}_{3}}) = ({{f}_{1}},{{f}_{2}},{{f}_{3}})$, то следует построить проекцию ${{P}_{3}}({{f}_{2}},{{f}_{3}},{{f}_{4}})$ трехмерного вектора (${{f}_{2}},{{f}_{3}},{{f}_{4}}$) на конус ${{V}_{3}}$, определенный в (4), и проверить включение $({{f}_{1}},{{P}_{3}}({{f}_{2}},{{f}_{3}},{{f}_{4}})) \in {{V}_{{3 + 1}}}$. Если это включение выполнено, то $({{f}_{1}},{{P}_{3}}({{f}_{2}},{{f}_{3}},{{f}_{4}}))$ – искомая проекция, так как выполнены условия леммы при $A = 0$, ${\mathbf{c}} = 0$ для ${{V}_{{1 + 3}}} = {{V}_{{3 + 1}}}$.

2. Если ${{P}_{3}}({{f}_{1}},{{f}_{2}},{{f}_{3}}) \ne ({{f}_{1}},{{f}_{2}},{{f}_{3}})$, то проекция ${{P}_{{3 + 1}}}({{f}_{1}},{{f}_{2}},{{f}_{3}},{{f}_{4}})$ есть результат аппроксимации фрагмента (${{f}_{1}}$, ${{f}_{2}}$, ${{f}_{3}}$, ${{f}_{4}}$) временного ряда фрагментом, значения которого есть линейная функция времени, в соответствии с условиями леммы при $B = 0$, $C = 0$, ${\mathbf{d}} = 0$, ${\mathbf{g}} = 0$.

Итак, если заданы параметры формы ${{V}_{f}}({{i}_{1}}, \ldots ,{{i}_{{2m}}})$ временного ряда, то построить наилучшее приближение любого временного ряда рядами заданной формы можно, основываясь на приведенных здесь результатах. Однако на практике интервалы выпуклости вверх и вниз могут быть непостоянными и меняться в известных пределах: $({{i}_{1}}, \ldots ,{{i}_{{2m}}}) \in I$, $I = \{ ({{i}_{1}}, \ldots ,{{i}_{{2m}}}):{{\underline i }_{k}} \leqslant {{i}_{k}} \leqslant {{\bar {i}}_{k}},\;k = 1, \ldots ,2m\} $. Интерес представляет выбор значений этих параметров, при которых анализируемый временной ряд наиболее близок к соответствующей форме ${{V}_{f}}({{i}_{1}}, \ldots ,{{i}_{{2m}}})$.

Задачу наилучшего приближения временного ряда сигналами формы ${{V}_{f}}({{i}_{1}}, \ldots ,{{i}_{{2m}}})$ выбором параметров формы назовем задачей морфологической фильтрации временного ряда. Формально это есть задача на поиск минимума

где I – множество допустимых значений параметров формы ${{V}_{f}}({{i}_{1}}, \ldots ,{{i}_{{2m}}})$.

4. ЗАДАЧА ФИЛЬТРАЦИИ СИГНАЛА ЗАДАННОЙ ФОРМЫ

Рассмотрим пример анализа временного ряда, при котором используется морфологическая фильтрация, описанная в разд. 2 настоящей статьи. На фиг. 1 точками отображена временная зависимость концентрации СО₂, измеренная с шагом 0.5 ч. По вертикали отложена концентрация в миллионных долях (ppm), по горизонтали – время. Видно, что фрагменты временного ряда длительностью примерно 24 ч схожи по форме, однако имеются вариации как в амплитуде, так и в длительности фрагментов. В целом, 24-часовые фрагменты хорошо описываются сигналами, математическая модель формы которых описана в предыдущем разделе настоящей статьи, однако имеются и отличия, которые тоже могут представлять интерес для исследователя.

Для описания суточного хода концентрации СО₂ была применена морфологическая фильтрация сигнала. Морфологическая фильтрация осуществляется следующим образом. Задаются ограничения на параметры формы, затем для каждого возможного значения параметров ${{i}_{1}}$ и ${{i}_{2}}$ строится проекция фрагмента временного ряда на множество фрагментов, выпуклых вверх, и выбираются такие их значения, которые минимизируют отличие анализируемого фрагмента от заданной формы. Затем для найденного значения ${{i}_{2}}$ и для каждого возможного значения ${{i}_{3}}$ строится проекция фрагмента временного ряда на множество фрагментов, выпуклых вниз. Далее проверяется условие сшивания фрагментов, выпуклых вниз и вверх. Условиe сшивания задается неравенствами (3). Если они не выполнены, корректируются значение найденной проекции и значение ${{i}_{2}}$ положения точки перегиба, далее эта процедура повторяется для каждого участка заданного типа выпуклости.

На фиг. 1 приведен результат морфологической фильтрации (сплошная линия). Хорошо видны повторяющиеся участки длительностью примерно 24 ч, отвечающие суточной изменчивости концентрации СО₂. Средний период, вычисленный по исследованному фрагменту временного ряда, составил 24.0 ч.

В результате применения алгоритма морфологической фильтрации к временному ряду концентрации СО₂ наряду с выделением суточной циклической составляющей из исходного временного ряда проведено исследование остаточной составляющей, представляющей собой разность исходного ряда и его аппроксимации, полученной в результате морфологической фильтрации. График значений остатков показан на фиг. 2; на нем явно угадываются повторяющиеся парные пики. Они были выделены из остаточной составляющей повторной морфологической фильтрацией с периодом, меньшим, чем 24 ч. Результат повторной морфологической фильтрации приведен на фиг. 2 в виде кривой оранжевого цвета.

Проверка остаточной составляющей на стационарность с помощью расширенного теста Дики–Фуллера [12] показала, что, в отличие от исходного ряда, ряд остатков удовлетворяет условию стационарности, следовательно, для дальнейшего выделения циклических компонент ряда с периодами в пределах суток возможно применение фурье-анализа.

Приближенное представление ряда остатков ${{X}_{t}}$ в виде частичной суммы ряда Фурье

где $t = 1,\;2,\; \ldots ,\;N$, $N$ – число элементов ряда, ${{a}_{i}}$, ${{b}_{i}}$ – коэффициенты Фурье, ${{T}_{i}}$ – период i-й гармоники, $i = 1,\; \ldots ,\;I$, I – число гармоник, позволяет обнаружить колебания ряда ${{X}_{t}}$ остаточной составляющей с периодами в пределах суток, вносящие наибольший вклад во временной ряд ${{X}_{t}}$, и оценить вклад i-й гармоники с помощью периодограммы – значений ${\text{Pe}}{{{\text{r}}}_{i}} = \frac{N}{2}({{a}_{i}}^{2} + {{b}_{i}}^{2})$, $i = 1,\; \ldots ,\;I$, а ее относительный вклад в суммарную периодограмму ряда $({{X}_{t}} - {{a}_{0}})$ оценить значением

В результате применения фурье-анализа к ряду остаточной составляющей показано, что наиболее существенный вклад вносят гармонические компоненты разложения Фурье (5) с периодами 11.9, 7.95 и 6 ч.

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

5. СТАТИСТИЧЕСКИЙ АНАЛИЗ ВРЕМЕННОГО РЯДА КОНЦЕНТРАЦИИ УГЛЕКИСЛОГО ГАЗА

Статистический анализ временных рядов концентрации углекислого газа также проводится на основе данных эколого-климатической станции “AsiaFlux” во Вьетнаме с 2011 по 2018 г. На фиг. 3 графически представлен временной ряд концентрации CO₂ за период с 2011 по 2018 г. На основе визуального анализа графика концентрации CO₂ выдвинем предположение о том, что представленный ряд содержит регулярную и случайную компоненты. Такое предположение подтверждается также заключениями метеорологов о механизмах формирования уровня концентрации CO₂ в атмосфере под воздействием различных природных и антропогенных факторов. При этом регулярная компонента включает основную тенденцию ряда динамики показателя (тренд), а также его сезонные и циклические составляющие, характеризующиеся различными периодами колебаний концентрации CO₂.

Визуальный анализ рядов динамики концентрации CO₂ позволяет также высказать предположение о том, что значительных изменений амплитуды колебаний концентрации CO₂ с 2011 г. по 2018 г. как в сторону заметного увеличения, так и в сторону уменьшения, не наблюдается, поэтому в качестве модели декомпозиции изучаемых рядов динамики предлагается аддитивная модель [13]:

где ${{S}_{t}}$ – постоянная сезонная составляющая, ${{C}_{t}}$ – циклическая составляющая, ${{u}_{t}}$ – тренд, определяющий основную тенденцию временного ряда, ${{\varepsilon }_{t}}$ – нерегулярная составляющая, $t = 1,\;2,\; \ldots ,\;N$, где $N$ – число элементов ряда.

При изучении основной тенденции изменчивости концентрации CO₂ за рассматриваемый временной период целесообразно на начальном этапе, предшествующем декомпозиции, провести процедуру сглаживания исходного ряда с целью повышения его устойчивости по отношению к выбросам различной природы. В результате сглаживания методом простой скользящей средней исходного ряда динамики концентрации CO₂${{X}_{t}}$, $t = 1,\; \ldots ,\;N$ (6), получаем ряд скользящих средних ${{\hat {X}}_{i}}$:

где K – интервал сглаживания, равный периоду колебаний концентрации CO₂, который будем задавать различным образом, учитывая периоды колебаний, вносящие наибольший вклад в ряд динамики концентрации CO₂. Если провести усреднение ежесекундных значений концентрации CO₂ не за полчаса, а за сутки, т.е. в качестве ${{X}_{t}}$, $t = 1,\; \ldots ,\;N$, рассмотреть среднесуточные значения уровней ряда, а интервал сглаживания задать равным 365 дней, в результате сглаживания получится ряд ${{\hat {X}}_{t}}$, $t = 1,\; \ldots ,\;N$, не содержащий влияния сезонных изменений. Тогда для определения постоянной сезонной компоненты ${{S}_{t}}$ следует усреднить разность значений ${{X}_{t}}$ и ${{\hat {X}}_{i}}$ для каждого дня $t = k + Kj$, $t = 1,\;2,\; \ldots ,\;N$, на протяжении восьмилетнего периода наблюдений за концентрацией CO₂, где $k$ – день года, j – год, k = 1, …, K, $j = 0,\; \ldots ,\;J - 1$: где $J = 8$ лет, K = 365 дней. Для дальнейшего изучения основной тенденции ряда (6), не зависящей от сезонных колебаний, следует вычесть из исходного ряда ${{X}_{t}}$, $t = 1,\; \ldots ,\;N$, ряд сезонной компоненты ${{S}_{t}}$, $t = 1,\; \ldots ,\;N$, и снова провести усреднение ряда ${{X}_{t}} - {{S}_{t}}$, $t = 1,\; \ldots ,\;N$, методом скользящих средних (7). Поскольку в этом случае характер основной тенденции может зависеть от выбора периода усреднения, проведем анализ циклических колебаний ряда динамики концентрации CO₂.

Наряду с временными рядами концентрации CO₂ на различных высотах над поверхностью Земли, важное значение для анализа климатических изменений имеют ряды динамики экстремальных (минимальных и максимальных) суточных значений концентрации CO₂ (фиг. 4). Проверка временных рядов концентрации CO₂ на стационарность с помощью расширенного теста Дики–Фуллера [12] показала, что в то время как ряды динамики нельзя считать стационарными, временные ряды экстремальных суточных значений концентрации CO₂ являются стационарными, поэтому для анализа их циклических колебаний применен фурье-анализ. В результате применения фурье-анализа к ряду минимальных суточных значений концентрации CO₂ выделены следующие колебания концентрации CO₂, вносящие наибольший вклад во временной ряд ${{X}_{t}}$: 11.37 мес, 3.73 г, 15.1, 6.5, 9, 7.6 и 3.5 мес (перечислены в порядке убывания значений периодограммы). Графическое представление циклической составляющей C_t временного ряда (6) показано на фиг. 5. При анализе максимальной суточной концентрации были выделены совпадающие по значениям, но отличающиеся по величине периодограммы периоды, приведенные далее также в порядке убывания значений периодограммы: 6.5, 15.1, 11.37, 7.6, 5.7 мес. Проверка остаточных составляющих ${{\varepsilon }_{t}}$ (5) рядов максимальных и минимальных суточных концентраций с помощью критерия Колмогорова–Смирнова показала, что ряды остатков подчиняются нормальному закону распределения.

Таким образом, с помощью морфологического анализа исходного временного ряда концентрации CO₂ и остаточной составляющей, представляющей собой разность исходного ряда и его аппроксимации, полученной в результате морфологической фильтрации, были выявлены суточные колебания концентрации CO₂, а также циклические составляющие с периодами в пределах суток: 12, 8 и 6 ч. Для определения циклических компонент, имеющих более длительные периоды, был применен фурье-анализ стационарных рядов минимальных и максимальных суточных значений концентрации CO₂. В результате фурье-анализа показано, что наибольший вклад в исследуемые за период с 2011 по 2018 г. временные ряды вносят гармоники около 1 г., 3.7 года, 6.5 и 9 мес.

6. ВЕЙВЛЕТ-АНАЛИЗ ВРЕМЕННОГО РЯДА КОНЦЕНТРАЦИИ УГЛЕКИСЛОГО ГАЗА

В случае, когда ряд не является стационарным и преобразовать его в стационарный не удается, оправдано применение вейвлет-анализа. Вейвлеты – это семейство функций, которые получаются из одной функции посредством ее сдвигов и растяжений по оси времени. С помощью вейвлет-преобразования функция рассматривается в виде разложения на колебания, локализованные по времени и частоте. Таким образом, вейвлет-анализ применяется для обработки нестационарных сигналов, так как отвечает специфике временных рядов, которые демонстрируют эволюцию во времени своих основных характеристик – среднего значения, дисперсии, периодов, амплитуду и т.д. [14] и [15].

Исследуемый временной ряд показателей концентрации СО₂ за полтора года, усредненных с периодом в 30 мин, в ходе работы был проверен на стационарность с помощью расширенного теста Дикки–Фуллера. Результат оказался предсказуем – ряды нестационарны, тем самым обосновано применение вейвлет-анализа. Временные ряды минимумов и максимумов на различных высотах, напротив, оказались стационарными, что, вообще говоря, не исключает возможность использования вейвлет-преобразования в качестве одной из методик исследования, в связи с этим вышеупомянутые ряды динамики также использовались в качестве анализируемых временных рядов.

Сущность непрерывного вейвлет-преобразования заключается в следующем. С помощью подходящего материнского вейвлета $\Psi (t)$ вычисляются вейвлет-функции ${{\Psi }_{{a,b}}}(t)$:

где параметр $a$, называемый параметром масштаба вейвлет-преобразования, принимает строго положительные значения и отвечает за ширину вейвлета; величина $b$ есть параметр сдвига, который определяет положение вейвлета на оси $t$.

Далее определяются вейвлет-коэффициенты $W(a,b)$:

где f(t) – анализируемый ряд динамики, * обозначает комплексное сопряжение [5].

После этого проводится качественный анализ картины вейвлет-коэффициентов и построение интегрального спектра:

который показывает наличие циклов в исходном временном ряде. Масштаб $a$ связан с координатами оси времени $t$ как $t = \frac{a}{{{{F}_{c}}}}$, где ${{F}_{c}}$ – центральная частота вейвлета.

Для исследуемых временных рядов в качестве материнских вейвлетов были выбраны вейвлет Морле и вейвлет Гаусса 7, дающие наиболее информативные и наглядные результаты [16]. Центральная частота данных материнских вейвлетов – 0.8125 и 0.6 соответственно (фиг. 6). Применение иных материнских вейвлетов давали результаты, принципиально не отличающиеся от приведенных в настоящей статье, но с менее выраженными частотно-временными особенностями.

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

Интегральный спектр, в отличие от картины вейвлет-коэффициентов, позволяет рассчитать величину периода цикла (фиг. 7б и 8б). Локальные экстремумы графика соответствуют максимальным коэффициентам на данном масштабе, показывая значения периода циклов.

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

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

Соответствующий интегральный спектр (фиг. 7б) также является достаточно размытым. Можно наблюдать пики в 23, 45, около 80 дней для минимума концентрации, а также общие и для показателей минимумов и максимумов пики в 1.04 и 1.5 года.

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

При анализе данных концентрации углекислого газа за полтора года с усреднением в полчаса также выявлено множество непродолжительных (вносящих вклад в ряд динамики на протяжении только части области исследования) цикличностей (фиг. 8а). В связи с чем для получения более информативных результатов необходимо анализировать отдельные участки временного ряда, которые следует выделить в согласовании с климатическими условиями местности.

7. ЗАКЛЮЧЕНИЕ

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

Метод морфологической фильтрации дает возможность преобразовать нестационарный ряд в стационарный, что позволяет провести дальнейшее исследование ряда методами фурье-анализа и статистики. В случае нестационарных рядов выделить циклы также позволяет применение вейвлет-анализа. В нашем случае – это 24, 12, 8 и 6 ч для периодов, не превосходящих суточные. Метод вейвлет-анализа позволил локализовать циклические участки с периодом в несколько дней: 23, 45, около 80 дней для минимума концентрации, а также общие и для показателей минимумов и максимумов на всех высотах пики в 1.04 и 1.5 года.

Периодичность в 1.04 года соответствует в целом годовому циклу изменения концентрации СО₂ и обусловлен сезонной динамикой осадков и фенологическими изменениями в развитии растительности, а в 24 ч соответствует суточному циклу. Объяснение природы циклов с другими периодами требует дальнейших исследований.

Таким образом, предложен набор методов, позволяющих выявить цикличность рядов в широком временном диапазоне.