Исследование природных ритмов хвойных деревьев

Актуальность проблемы

С середины XIX века средняя глобальная температура Земли увеличилась на 0,75 ± 0,2 °C градуса (IPCC, 2007).

Причина этого явления во многом до сих пор не ясна. Часть учёных полагает, что сегодняшнее потепление — это закономерное явление и оно связано с естественной цикличностью процессов. Поскольку климатические изменения подобного масштаба наблюдались и ранее, не вызывая серьезных изменений на уровне экосистем, нынешнее потепление можно рассматривать с точки зрения естественной изменчивости климата Земли.

Противники этой теории считают, что глобальное потепление — это уникальное явление, которое может принести непредсказуемые последствия и серьёзно изменить климат. Они утверждают, что потепление происходит в связи с существенным увеличением концентрации парниковых газов в атмосфере (углекислота, метан, закись углерода), которое, в свою очередь, вызвано человеческой деятельностью (антропогенный эффект). Предполагается, что эти повышение температуры будет не только превышать естественные природные колебания, но и происходить в 15−40 раз быстрее, чем аналогичные изменения в прошлом (IPCC, 2007).

Согласно разным сценариям выбросов парниковых газов, удвоение концентрации эквивалента углекислоты в атмосфере по сравнению с доиндустриальным уровнем произойдет между 2025—2050 гг., что изменит глобальную среднюю температуру воздуха от 1,5 до 4,5°С. Средняя скорость увеличения глобальной средней температуры в следующем столетии может составить примерно 0,3°С за десятилетие. Поверхность суши нагреется быстрее по сравнению с поверхностью океана.

Ожидается, что это потепление будет неравномерным: в тропических районах оно составит половину глобального среднего, а наибольшим будет в полярных районах, где превысит среднюю глобальную величину примерно в 2 раза и в основном произойдет в зимние месяцы (IPCC, 2007).

Рис. 1. Тренд температурной изменчивости поверхности Земли в градусах Цельсия за десятилетие: вверху — начиная с 1870 г., внизу — с 1950 г (Jones et. al., 2001)

Возможности прогноза естественной и антропогенной изменчивости климата, часто ограничивается в связи с короткими рядами метеорологических данных. Как правило, их длина не превышает 30−50 лет, за редкими исключениями — 100 лет. Косвенные индикаторы климата — ежегодный прирост древесины (ширина колец), толщина снежного покрова на ледниках, кораллы и прочие — могут в какой-то мере заменить инструментальные наблюдения в прошлом и продолжить ряды. Дендрохронологический метод в этом отношении является одним из самых перспективных: эти данные имеют годовое и сезонное разрешение, поддаются калибровке с помощью инструментальных метеорологических данных и потенциально доступны во всех районах, где растут деревья.

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

Цель работы

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

Для достижения цели необходимо:

1. собрать и структурировать данные по приросту древесины на севере Европейской части России и построить сводные хронологии ширины и плотности древесины

2. провести корреляционный анализ между площадками для объединения их в региональные хронологии

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

4. с помощью вейвлет-анализа выявить периодичность в хронологиях и сравнить эти результаты с результатами, полученными методами преобразований Фурье и СВАН-диаграмм

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

6. оценить пространственную и временную изменчивость природных ритмов за последние 300 лет на субрегиональном уровне

Стратегия исследования

На схеме ниже схематично показана последовательность действий, которой мы придерживались в этой работе.

Хронологии, создаваемые до 1990х гг. на европейской части России и на прилегающих территориях за немногими исключениями (Шиятов, 1986), не подвергались хронологическому контролю (перекрестное датирование) (Fritts, 1976, Шиятов и др., 2002) и поэтому не могут использоваться для дальнейшего анализа. Хронологии, создаваемые на основе археологической древесины, не сведены в региональные, недоступны широкому пользователю, и используются лишь для датирования археологических и архитектурных объектов (Черных, 1989).

В 1990;х гг. благодаря исследованиям международной группы специалистов (Schweingruber, Briffa, 1996; Ваганов и др., 1996; Briffa et al., 2001) удалось получить большой массив дендрохронологической информации по ширине и плотности колец на северной границе леса в Евразии, эти данные доступны и хранятся в международном банке данных (www.cdc.noaa.gov). Для юга Русской равнины таких данных нет, поскольку в этом районе крайне трудно найти высоковозрастные деревья и погребенную древесину. В окружающих Русскую равнину горах построены древесно-кольцевые хронологии для Полярного и Приполярного Урала (Шиятов, 1986) и Кольского полуострова (Hiller et al., 2001, Шумилов и др., 2007).

Многочисленные древесно-кольцевые хронологии лиственницы (Larix sibirica, L. dahurica) и ели (Picea abies, P. obovata) вдоль северной границы леса между 42°E и 152°E, за последние 400−600 лет показывают хорошую взаимную согласованность (Schweingruber, Briffa, 1996; Vaganov et al., 1996; Briffa et al., 2001; 2004), поэтому они были сведены в несколько региональных хронологий. В частности, одна из них представляет европейский сектор (хронология NEUR) и включает (помимо российских) скандинавские хронологии (Briffa et al., 2001). Как ширина кольца, так и максимальная плотность древесины этих хронологий положительно коррелируют с температурой теплого периода (апрель-сентябрь r~0.7). Однако, эти реконстркукции были ограничены 1960;м годом, поскольку для более позднего периода был обнаружен эффект рассогласования дендрохронологических и климатических показателей.

Для Урала (Шиятов, 1986) и Кольского полуострова (Hiller et al., 2001) реконструкции температур теплого периода охватывают гораздо более длительный период, поскольку хронологии включают погребенную древесину. Более того, для реконструкции долгопериодных трендов в обоих районах используются также данные по динамике верхней границы леса. Работы на Кольском полуострове позволили создать несколько длительных хронологий (Hiller et al., 2001). Хронология показывает значимую положительную корреляцию с солнечной активностью и лишь невысокие положительные коэффициенты (до 0,25) с температурами с августа по ноябрь. Сходые результаты были получены и Распоповым с соавторами для сосны для этого же района (Raspopov et al. (2004).

Наиболее систематические и исчерпывающие по интересующему нас вопросу были получены К. Бриффой с соавторами (Briffa et al., 2004), Эти авторы хронологии плотности древесины осредняли по сетке 5×5 градусов и коррелировали со среднемесячной температурой осредненной для того же квадрата. Функции отклика для каждой локальной хронологии были рассчитаны с опорой на данные о среднемесячных температурах (Jones et al., 1999) и осадках (New et al., 2000) по регулярной сетке, для предшествующего формированию кольца и текущего года (за период 1881—1984 гг.). Анализ показал (рис.2) (Briffa et al., 2004), что наиболее тесные корреляции в целом обнаруживаются для максимальной плотности годичных колец и температур в период с апреля по сентябрь.

Рис. 2. Функция отклика для максимальной плотности древесины хвойных на северной границе леса (Briffa et al., 2004)

Реконструкции осадков для северных районов по годичным кольцам деревьев крайне редки, поскольку лимитирующим фактором в этих условиях является в подавляющем большинстве случаев летняя температура. Оригинальную попытку реконструировать осадки на основе дендрохронологических данных Швейнгрубера на Полярном Урале предпринял Ю. М. Кононов с соавторами (Kononov, Ananicheva, Willis, 2005). Они в частности, показали, что ширина колец ранней древесины положительно коррелирует с осадками холодного периода, однако насколько устойчива эта связь пока остается неясным из-за ограниченности привлеченного Кононова с соавторами материала и неясного экологического обоснования такой связи.

Таким образом, несмотря на подробный анализ данных дендрохронологии на севере Евразии, проведенный нашими предшественниками (Ваганов и др, Briffa et al., 2004), многие вопросы, особенно касающиеся региональной климатической изменчивости в этом районе, остались неисследованными. В частности, неизвестно насколько отличалась динамика ширины и плотности колец с запада на восток — на Кольском полуострове, в центральной части севера Русской равнины и на Урале, каковы основные ритмы этой изменчивости и насколько они устойчивы во времени за последние 300 лет, действительно ли по ширине ранней древесины возможно реконструировать осадки холодного периода в этом регионе и др. Эти проблемы будут, в частности, рассмотрены в этом исследовании.

Глава 1. Основы дендрохронологии

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

Второе обстоятельство состоит в том, что способность дендрохронологических реконструкций воспроизводить долгопериодные колебания климата ограничена в связи с тем, что при создании сводных хронологий приходится проводить индексирование серий, необходимое для вычленения возрастного тренда прироста (Шиятов и др., 2000). Это, однако, не мешает впоследствии оценить на основе древесно-кольцевых хронологий колебания климатических параметров межгодовой и внутривековой продолжительности, в частности, частоту экстремумов, периодичность колебаний, тренды и некоторые другие статистические характеристики.

Методы дендрохронологии в реконструкции климата

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

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

Основные принципы дендрохронологии

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

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

Принцип чувствительности: в благоприятных для роста деревьев районах и местообитаниях формируются широкие годичные кольца. При этом у таких деревьев хорошо выражены изменения прироста с возрастом, а величина прироста между соседними годами колеблется в незначительных пределах. Подобная последовательность называется «благодушной». В неблагоприятных для произрастания деревьев условиях кольца прироста узкие, их ширина значительно колеблется от года к году, возрастная кривая роста выражена слабо. Часто наблюдается выпадение колец. Такие серии колец называются «чувствительными». Они свидетельствуют о том, что на рост деревьев большое влияние оказывают факторы внешней среды.

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

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

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

Достоинства и недостатки метода

Как и любой метод, дендрохронология обладает рядом достоинств и недостатков.

Достоинства дендрохронологических методов:

· Высокая разрешающая способность древесно-кольцевых хронологий (год и сезон года).

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

· Возможность получения длительных и непрерывных хронологий (сотни и тысячи лет).

· Возможность выявлять колебания различной длительности в изменчивости различных характеристик прироста деревьев и определяющих прирост факторов внешней среды.

· Возможность получения информации для огромных территорий суши, где растут или росли в прошлом древесные растения.

Недостатки дендрохронологических методов:

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

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

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

· Сходная реакция прироста древесины на различные факторы среды.

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

Измерение характеристик древесины

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

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

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

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

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

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

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

Рис. 3. Пример построения древесно-кольцевой хронологии длительностью 1950 лет на основе использования образцов древесины, взятых с живого дерева и двух деревьев, отмерших в 1228 и 897 годах н.э. соответственно.

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

Дендроклиматические материалы и методы

В этой работе использовались результаты измерений ширины ранней древесины и плотности колец отдельных образцов из международного банка дендрохронологических данных (www.cdc.noaa.gov).

Ширина ранней древесины

Плотность древесины

Рис. 4. Наглядное представление ранней ширины и плотности древесины

Процедура дендроклиматического анализа состоит в следующем. Образцы древесины (керны и спилы) отбираются на некоторой ограниченной территории с однородными орографическими, почвенными и микроклиматическими условиями. Деревья на площадках выбираются наиболее старые, отдельно стоящие, не угнетённые, не повреждённые, так как основной фактор, влияющий на прирост таких деревьев, может быть неклиматическим (Schweingruber, 1988; Cook, 1990). Для уменьшения вероятности выпавших или ложных колец в кернах, с одного ствола отбираются обычно по два керна по двум произвольным радиусам на высоте 1−1,5 м от поверхности земли.

Дальнейшая работа с образцами древесины проводится в лабораторных условиях. Керны наклеиваются на специальную деревянную основу, после чего керн полируют для увеличения изображения контрастности колец (рис.5). Затем, зная точную календарную дату отбора образца, проводится разметка и предварительная датировка колец (Шиятов и др., 2000).

Рис. 5. Образцы, подготовленные для дендрохронологического анализа

С помощью полуавтоматического комплекса состоящего из микроскопа и подвижного измерительного стола, измерялась ширина годичных колец с точностью до ±0,01 мм. Величина ширины колец автоматически фиксируется компьютерной программой TSAP (Time Series & Analysis Presentation — программа анализа и перекрестной датировки дендрохронологических серий и графического представления результатов обработки данных), которая помогает решить широкий спектр задач, связанных с анализом древесных кернов и спилов.

Метод перекрёстного датирования позволяет с высокой точностью определить календарную дату выпадающих или ложных колец. Кроме того, вероятность выпадения колец увеличивается с возрастом. Программа COFECHA позволяет произвести контроль качества полученных древесно-кольцевых серий. Для проверки правильности датировки программа COFECHA преобразовывает серии измерений с помощью короткого кубического сплайна, а затем сравнивает все полученные серии и выявляет проблематичные участки (Holmes, 1983). Каждая индивидуальная серия делится на сегменты по 50 лет, каждый сегмент имеет общий интервал продолжительностью 25 лет с соседним сегментом и т. д. Значения коэффициентов корреляции определяются программой для каждого сегмента серии, что наиболее точно может указать на ошибку датирования. Статистически значимым считается коэффициент корреляции, значение которого более, чем 0,3281. Кроме того, это программа позволяет определить абсолютный возраст погребённой древесины и образцов сухих деревьв относительно имеющейся хронологии.

Рис. 6.Значения коэффициентов корреляции некоторых кернов со всей хронологией Кавказа

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

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

Рис. 7. Пример индексирования — избавления от возрастного тренда

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

Глава 2. Методы обработки нестационарных сигналов

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

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

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

Существует несколько подходов для анализа нестационарных сигналов такого типа. Первый — локальное преобразование Фурье (Short-Time Fourier Transform), оно же называется оконным преобразованием. В этом случае, мы работаем с нестационарным сигналом, как со стационарным, предварительно разбив его на сегменты (окна), статистика которых не меняется со временем (в данной работе это преобразование будем называть СВАН — спектрально-временной анализ), но у этого метода также есть ряд небольших недостатков, которые мы обсудим позже. Второй подход — вейвлет преобразование. В этом случае нестационарный сигнал анализируется путем разложения по базисным функциям, полученным из некоторого прототипа путем сжатий, растяжений и сдвигов. Функция-прототип называется материнским, или анализирующим вейвлетом.

Краткий обзор преобразования Фурье

Традиционно для анализа временных рядов используется преобразование Фурье — преобразование функции, превращающее её в совокупность частотных составляющих. Более точно, преобразование Фурье — это интегральное преобразование, которое раскладывает исходную функцию по базисным функциям, в качестве которых выступают синусоидальные (или мнимые экспоненты) функции, то есть представляет исходную функцию в виде интеграла синусоид (мнимых экспонент) различной частоты, амплитуды и фазы. В более общей форме записи:

где

— амплитуды гармонических колебаний соответствующей частоты, которые определяются по формуле:

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

Преобразования:

· не отличают сигнал, являющийся суммой двух синусоид, от ситуации последовательного включения синусоид

· не дают информации о преимущественном распределении частот во времени

· могут дать неверные результаты для сигналов с участками резкого изменения.

Исследуемые ряды также далеко не всегда удовлетворяют требованию периодичности и более того, как правило, заданы на ограниченном отрезке времени.

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

Одним из таких методов является обобщенный метод Фурье (локальное (оконное) преобразование Фурье). Этот метод состоит из следующих этапов:

1. в исследуемой функции создается «окно» — временной интервал, для которого функция f (x) 0, и f (x)=0 для остальных значений,

2. для этого «окна» вычисляется преобразование Фурье,

3. «окно» сдвигается, и для него также вычисляется преобразование Фурье Пройдя таким «окном» вдоль всего сигнала, получается некоторая трехмерная функция, зависящая от положения «окна» и частоты.

Но данный подход позволяет определить факт присутствия в сигнале любой частоты, и интервал ее присутствия. Это значительно расширяет возможности метода по сравнению с классическим преобразованием Фурье, но существуют и определенные недостатки. Согласно следствиям принципа неопределенности Гейзенберга в данном случае нельзя утверждать факт наличия частоты w₀ в сигнале в момент времени t₀ — можно лишь, что спектр частот (w₁, w₂) присутствует в интервале (t₁, t₂). Причем разрешение по частоте (по времени) остается постоянным вне зависимости от области частот (времен), в которых производится исследование. Поэтому, если, например, в сигнале существенна только высокочастотная составляющая, то увеличить разрешение можно только изменив параметры метода.

В данной работе такой метод мы будем называть СВАН-преобразованием.

В качестве метода, не обладающего подобного рода недостатками, был предложен аппарат вейвлет анализа.

Основные положения вейвлет-преобразования

Альтернатива преобразованию Фурье для исследования временных (пространственных) рядов с выраженной неоднородностью стал метод, разработанный в 80-х годах под названием вейвлет-анализ.

Вейвлеты (от англ. wavelet), всплески — это математические функции, позволяющие анализировать различные частотные компоненты данных.

Вейвлеты — это семейство функций, которые локальны во времени и по частоте («маленькие»), и в которых все функции получаются из одной посредством её сдвигов и растяжений по оси времени (так что они «идут друг за другом»).

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

Вейвлет-преобразования обычно делят на дискретное вейвлет-преобразование (ДВП) и непрерывное вейвлет-преобразование (НВП). ДВП обычно используется для кодирования сигналов, в то время как НВП для анализа сигналов. В результате, ДВП широко применяется в инженерном деле и компьютерных науках, а НВП в научных исследованиях.

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

Обычно, функция-вейвлет обозначается буквой ш.

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

При конструировании базисной анализирующей функции должны выполняться следующие необходимые условия:

· Локализация — вейвлет должен быть локализован вблизи нуля аргумента как во временном, так и в частотном пространстве. Нулевое среднее: .

· Вейвлет должен быть знакопеременной функцией. Ограниченность: .

· Вейвлет должен достаточно быстро убывать по временной (пространственной) переменной.

Наибольшей популярностью пользуются вейвлеты:

Непрерывное вейвлет-преобразование (НВП) строится с помощью непрерывных масштабных преобразований и переносов вейвлета с произвольными значениями масштабного коэффициента a и параметра сдвига b:

(1)

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

Вейвлет-преобразование обратимо для функций f из L²®

.

Итак, у нас имеется некоторая функция f (t), зависящая от времени. Результатом ее вейвлет-преобразования будет некоторая функция W (a, b), которая зависит уже от двух переменных: от времени и от частоты (обратно пропорционально). Для каждой пары a и b рецепт вычисления вейвлет преобразования следующий:

1. Функция вейвлет растягивается в a раз по горизонтали и в 1/a раз по вертикали.

2. Далее он сдвигается в точку b. Полученный вейвлет обозначается .

3. Производится усреднение в окрестности точки a при помощи .

Спектр вейвлет-преобразования одномерного сигнала представляет поверхность в трехмерном пространстве. Обычно изображение спектра выполняется путем проектирования линий постоянного уровня поверхности на плоскость с переменными: параметрами сдвига (по оси абсцисс) и масштабом (по оси ординат), с градиентной заливкой оттенками серого цвета между линиями.

В результате «вырисовывается» наглядная картина, иллюстрирующая частотно-временные характеристики сигнала. По оси абсцисс откладывается время, по оси ординат — частота (иногда log (T), где T=1/a — период). А абсолютное значение вейвлет преобразования для конкретной пары a и b определяет цвет, которым данный результат будет отображен (чем в большей степени та или иная частота присутствует в сигнале в конкретный момент времени, тем более выраженный будет оттенок). На рисунке 8 показан пример вейвлет-анализа. Чётко видно, что на протяжении всего времени действия сигнала преобладают две частоты.

Рис. 8. Пример вейвлет-анализа. Внизу — исходный двухчастотный сигнал, вверху — вейвлет-анализ, показывающий присутствие сигналов с частотой примерно 0,8 и 3,5. По оси абсцисс время, по оси ординат — частота Таким образом, любая функция из LІ® может быть представлена суперпозицией масштабных преобразований и сдвигов базисного (материнского) вейвлета с коэффициентами, зависящими от масштаба (частоты) и параметра сдвига (времени).

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

Вейвлеты, используемые в данной работе:

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

Получив вейвлет-спектр, можно рассчитать полную энергию сигнала:

и глобальный спектр энергии — распределение полной энергии по масштабам частоты (скейлограмму вейвлет-преобразования) Скейлограмма соответствует спектру мощности Фурье-преобразования сигнала, сглаженному на каждом масштабе спектром Фурье анализирующего вейвлета:

где знак ^ обозначает Фурье-образ функции.

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

где

Дискретное вейвлет-преобразование

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

Базис нормирован, если вейвлет имеет единичную норму.

Вейвлет называется ортогональным, если семейство представляет ортонормированный базис функционального пространства L²®, т. е.. В этом случае любая функция может быть представлена в виде ряда:

Где Спектрально-временной анализ (СВАН) Как уже говорилось выше, в данной работе будет так же использоваться спектрально-временной анализ, он же — графическое представление оконного преобразования Фурье, то есть процедура спектрального анализа временных рядов в скользящем временном окне (предварительно исключается полиномиальный тренд, порядок которого соответствует той кратности, с которой заданное временное окно укладывается в полной длине ряда). Длина окна задается исходя из требуемой детальности и спектрального состава процесса (она должна быть кратна периодам исследуемых ритмов). Мы будем задавать длину 25% от длины ряда. Результат изображается в виде спектрально-временной диаграммы. На оси абсцисс откладывается календарное время (в годах), соответствующее центру скользящего временного окна. На оси ординат откладываются величины частот в циклах в единицу времени (год). Каждый вертикальный столбец представляет собой амплитудный спектр Фурье, рассчитанный в заданном скользящем временном окне. Более сильная зачерненность на диаграммах соответствует большим спектральным амплитудам. О величинах амплитуд судят по шкале уровней, помещенной около диаграммы. Повторяемость или устойчивость доминирующих ритмов выражается в виде протяженных зачерненных полос.

Но у этого анализа есть довольно существенный недостаток: в данном анализе чем шире выбираемое окно, тем больше «отрезается» значений от начального ряда после преобразования. Но, чем уже окно, тем хуже определяются частоты, доминирующие в данный момент времени. Например, если длина ряда 300 лет, а ширина окна 50 лет, то в поучившейся диаграмме будет отсутствовать по 25 лет с каждой стороны ряда.

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

Для оценки уровня упорядоченности процессов вычисляется т.н. «параметр хаотизации»:

где K = L/2+1 и L — длина скользящего временного окна анализа (количество отсчетов), W (, t) — амплитудный спектр на частоте и времени t. Нулевое значение параметра соответствует наличию одной спектральной линии (т.е. ряд состоит из одной синусоиды), а единичное — случаю белого шума. Таким образом, области повышенных значений параметра могут свидетельствовать о моментах повышенной «хаотизации» процесса (например, при распадной неустойчивости с образованием хаотических колебаний), а минимальные значения параметра — большей его «упорядоченности» (например, при синхронизации различных гармоник на единые предельные циклы).

Что касается достоверности определения исследуемых параметров, например, амплитуд и периодов циклических вариаций, зависит от погрешностей исходных временных рядов, дисперсии и статистического смещения рассчитываемых оценок параметров, разрешающей способности анализа. Иначе говоря, характеристика достоверности данных основывается на стандартных приемах доверительного оценивания. Как правило, мы стремимся определить спектрально-временную структуру процессов. Поэтому для нас наиболее важно знать, насколько достоверны выделяемые ритмические составляющие на СВАН-диаграммах. Интерпретируя их, мы оцениваем, во сколько раз амплитуды ритмических компонент превышают шум и сколько периодов укладывается в интервал времени, когда прослеживается та или иная выделяемая гармоника. Если амплитуда гармоники превышает амплитуду шума в несколько раз и выделяется достаточно протяженная полоса (в которую укладывается хотя бы несколько периодов), подтверждаемая наблюдениями в другие интервалы времени или на сходных объектах, то можно уверенно говорить о существовании в этом процессе данного ритма.

Сравнение методов анализа

Для того, чтобы понять, как все эти методы работают, наглядно представим результат обработки синусоидного сигнала каждым из них (рис.9)

Рис. 9. Сравнение методов анализа.

На рисунке 9 представлен синусоидный сигнал, частота которого увеличивается со временем, пусть для простоты единицы по оси x будут года. Как видно из рисунка, Фурье-преобразование даёт информацию о наличии той или иной частоты в данном сигнале, но при этом, когда именно эти частоты присутствовали нам не известно. В то же время СВАНи вейвлетанализы дают чёткую картину динамики изменения частотных характеристик во времени. На обеих рисунках видно, что доминирующий сигнал периодом примерно 120 лет, затем сменился на период длиной 35 лет, а в конце — примерно 5-ти летний период. Так же в СВАН-диаграмме по степени яркости всплесков видно, что амплитуда всех трёх составляющих сигналов одинакова. Но, как говорилось выше, в СВАН-диаграмме есть недостатки — потеря части ряда при преобразовании: это видно на рисунке 9 — в начале и в конце ряда отрезается несколько значений. Таким образом, можно сделать вывод, что вейвлет-преобразование и СВАН-анализ более информативны, по сравнению с преобразованием Фурье. Однако, для более детального анализа ряда всё же требуется Фурье-спектр, чтобы видеть, какие именно частоты доминируют в хронологии. Например, по данным Фурье-анализа видно, что две доминирующие частоты находятся в интервале от 0 до 0,05 [циклов/год]. В этом случае СВАН-анализ логичнее проводить в диапазоне именно этих частот (как и было сделано), чтобы более «крупно» видеть получившуюся картинку.

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

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

Глава 3. Ритмика природных процессов, отраженная в ширине и плотности колец на севере Русской равнины

Построение локальных хронологий и их выбор для включения в региональные; преобразование локальных хронологий в региональные

На рисунке 10 представлена северная часть европейской территории России и положение всех площадки, по которым в международном банке данных (www.cdc.noaa.gov) представлены первичные материалы по ширине и плотности колец. Площадки 1,3,5,6,8,22,27 сразу были отсеяны, так как по ним есть только данные по ширине колец, а данные по максимальной плотности отсутствуют.

Рис. 10. Положение дендрохронологических площадок на севере Русской равнины (www.cdc.noaa.gov)

Из набора измерений ширины и плотности колец отдельных образцов необходимо было составить сводные хронологии. Для этого мы использовали программу ARSTAN, которая исключает возрастной тренд путём индексирования ряда, а также исключает эффект эндогенности (внутренние причины) и помех со стороны других деревьев (Cook, Holmes, 1986). Пример такого индексирования представлен на рисунке 7. Таким образом, мы получили по 20 хронологий (всего 27 площадок минус 7 — те, по которым нет необходимых нам данных) максимальной плотности и ранней ширины древесины. После этого, с помощью коэффициента корреляции Пирсона мы рассчитывали взаимную корреляцию между площадками, чтобы выяснить, насколько сходны полученные хронологии между собой. Площадки, коэффициенты корреляции которых с остальными оказались незначимыми были исключены из дальнейшего анализа. По ширине колец были исключены площадки №№ 10,11,15,20, по максимальной плотности — №№ 10,13,15,20,24 (коэффициент корреляция ниже 0,3). В результате этой процедуры осталось 15 и 16 площадок, соответственно, которые затем Були объединены в несколько региональных хронологий.

Таким образом, по этому критерию площадки максимальной плотности разбились на 3 группы, из которых было составлено соответственно 3 региональных хронологии, а для ранней ширины — на 4. Объединение проведено путем погодичного осреднения индексов локальных хронологий. На рисунках 11 и 12 показано, как именно локальные площадки были объединены в региональные хронологии, под рисунком показаны коэффициенты корреляции и стрелками указано, к какой именно части территории относятся те или иные коэффициенты.

Коэффициент корреляции четырёх региональных хронологий по ширине равен от 0,27 до 0,56, а для трёх рядов по максимальной плотности — 0,37−0,67. Все коэффициенты корреляции значимы на 95% уровне. Это означает, что и ширина и плотность колец отражают некий общий климатический сигнал, причем каждый параметр — свой сигнал (ширина и плотность также скорелированы между собой, но с более низким коэфициентом корреляции).

Максимальная плотность несколько более тесно скоррелированый параметр и это дает основание предполагать, что ее связь с климатическими параметрами будет более заметной. Из работ наших предшественников (Ваганов и др., 1996, Briffa et al., 2004) известно, что максимальная плотность коррелирует с температурой теплого периода, а температура сама по себе — характеристика менее изменчивая в пространстве, чем осадки. Наш анализ косвенно подтвержадет этот вывод.

Таким образом у нас получилось всего 7 региональных хронологий (см. рис.13), которые далее исследуем с помощью вейвлет-, СВАНи Фурье-анализа. Выявим доминирующие частоты в данных рядах. Кроме того, изучим связь метео данных с нашими рядами.

Рис. 11. Разбиение площадок на региональные хронологии по ширине ранней древесины.

Региональные хронологии соответсвенно относятся к Кольскому полуострову, Карелии, северной части Русской равнины и Уралу. Внизу — соответствующие коэффициенты корреляции

Рис. 12. Разбиение площадок на региональные хронологии по максимальной плотности древесины.

Региональные хронологии соответственно относятся к Кольскому полуострову и Карелии, северной части Русской равнины и Уралу. Внизу — соответствующие коэффициенты корреляции

Рис. 13. хронологии ранней ширины и максимальной плотности колец хвойных на северной границе леса на Русской равнине. По оси ординат — индексы ширины и плотности колец

Корреляция хронологий с метеорологическими данными. Выявление доминирующего климатического сигнала

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

Рис. 14. Расположение метео станций использванных для анализа (синие точки)

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

Результаты корреляционного анализа представлены на рисунках 15,16,17. Согласно этим результатам, максимальная плотность древесины как на уровне локальных хронологий, так и для региональных, устойчиво коррелирует с температурой летнего периода (май-август) (рис.15). Коэффициент корреляции трёх региональных хронологий с данными летней температуры ближайших метеостанций и коэффициент отдельных локальных с этими же станциями составляет от 0,22 до 0,78) с уровнем значимости 95%. Таким образом, наши выводы согласуются с выводами К. Бриффы о корреляции максимальной плотности региональных хронологий для Сибири и севера Европы (включая Скандинавию) с температурой лета (см. рис.2).