автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.16, диссертация на тему:Новый подход к математическому моделированию временных рядов в экологии

На правах рукописи

Калуш Юрий Александрович

НОВЫЙ ПОДХОД К МАТЕМАТИЧЕСКОМУ МОДЕЛИРОВАНИЮ ВРЕМЕННЫХ РЯДОВ В ЭКОЛОГИИ

05.13.16 — применение вычислительной техники, математического моделирования и математических методов в научных исследованиях (в экологии)

Автореферат

диссертации на соискание учёноЗ степени кандидата физико-математических наук

Красноярск — 1998

Работа выполнена в лаборатории геоинформатики и моделирования процессов Тувинского института комплексного освоения природных ресурсов Сибирского отделения Российской академии наук

Научный руководитель: кандидат физико-математических наук, старший научный сотрудник ТувИКОПР СО РАН Логинов В.М.

Официальные оппоненты:

доктор физико-математических наук, профессор

Белолипецкий В.М.; кандидат физико-математических наук, профессор Джансеитов К.К.

Ведущая организация: Институт биофизики СО РАН (г. Красноярск)

Защита состоится 1 октября 1998 г. в 'часов на заседании диссертационного совета К 064.54.01 в Красноярском государственном техническом университете по адресу: г. Красноярск, ул. Киренского, 26

С диссертацией можно ознакомиться в библиотеке университета

Автореферат разослан 1 сентября 1998 г.

Учёный секретарь /И—'• диссертационного совета, ---дч Н.Г. Кузьменко

ОБЩАЯ ХАРАКТЕРИСТКА РАБОТЫ

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

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

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

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

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

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

числе BP метеорологического мониторинга в Туве и Северо-Западной Монголии (Улангом), рядов, связанных с глобальным изменением климата, ВР динамик численности популяций, всего более 50 рядов. Показана адекватность расчетных моделей фактическим данным.

Апробация работы. Результаты исследований докладывались на международном рабочем совещании "Комплексное изучение аридной зоны Центральной Азии" (Кызыл, 1994); Международном симпозиуме «Эксперимент Убсу-Нур» (Улангом, Монголия, 1995); Девятой Всесибирской школе по методам вычислительной математики (Шушенское, 1995); Семинаре по вопросам развития и применения геоинформационных технологий (г. Новосибирск, 1996); Втором сибирском конгрессе по прикладной и индустриальной математике (ИНПРИМ-96, Новосибирск, 1996); III межреспубликанском симпозиуме "Оптика атмосферы и океана" (Томск, 1996); III Международной научной конференции "Природные условия, история и культура Западной Монголии и сопредельных регионов" (Ховд, Монголия, 1997).

Публикации. Основные материалы диссертации опубликованы в семи научных работах [1—7]. Получено авторское свидетельство о государственной регистрации программы для ЭВМ [8].

Структура работы. Диссертация состоит из введения, трех глав, заключения, приложения и списка цитируемой литературы, насчитывающего 106 ссылок. Содержание диссертации в 101 страницу печатного текста иллюстрировано 38 рисунками.

СОДЕРЖАНИЕ РАБОТЫ

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

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

Рве. 1. Примеры временных радов.

1а) ВР "Глобальные изменения температуры". 16) ВР "Концентрация углекислого газа в атмосфере". [Доклад Гринпис "Глобальное потепление" — М.: Изд-во МГУ, 1993. —272 с.]

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

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

функции*^ происходит по "закону эффективного луча" /^(+)(/) на участках возрастания и Г(/) на участках убывания получаем, что ЪРх(1) на временах наблюдения может быть аппроксимирован некоторым эффективным временным рядом ^(+>(/) + -Проведенный в работе анализ распределения вероятностей значений углов раскрытия "веера" для участков монотонности показывает, что указанный разброс для реальных ВР, если произвести преобразование подобия, приводящее временные ряды к единому масштабу, весьма невелик и не превышает половины радиана. Приводится алгоритм выбора эффективного луча, как средневзвешенного от значений рассматриваемых лучей.

Рис. 2.

а. Схематический временной ряд. б,в. «Веера лучей», отображающие периоды роста и убывания ВР

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

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

на выделении составляющих временного ряда (таких как тренд, периодическая составляющая, шум), а также более современных методов моделирования ВР, не требующих указанного разбиения.

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

Во второй главе дается математическая формулировка задачи. Пусть вектор х(0 характеризует состояние рассматриваемой экосистемы (вектор состояния), и на временах / е Т.<+} компонентами х,(/) вектора состояния являются некоторые возрастающие функции, подчиняющиеся динамическому уравнению:

* = Г,{*'(я0),0' 0)

а на временах Г,'"' = Т \ Г,"' компоненты вектора состояния *,(/) описываются кусочно-у(5ывающими функциями, являющимися решением некоторого другого дифференциального уравнения:

х, = (2)

где Т— время наблюдения за поведением исследуемой экосистемы. Вводя индикаторную функцию а (О такую, что а (0=1, если * € Г,(+> И « (0 =-1, если / е7! , систему (1), (2) можно редуцировать к уравнению вида:

х, =7[/Г (Х.01 + -О (0[^<+) (х.О-^"' (*> 01,(3)

с начальным условием (/ = /„) = . Из уравнения (3) следует, что результирующее поведение системы, описываемое вектором состояния х, представляет собой результат перемешивания двух динамик. В зависимости от свойств множеств Г,(+) и Г.<_), на которых задается функция , это перемешивание носит детерминированный,

стохастический или смешанный характер. Привлекая изложенную выше качественную схему, примем, что в качестве функций F/*'(*(0, /) и F(("'(jr(0. ') служат производные по времени / от функций F'*' и Fполученных определенным усреднением изучаемого BP x,{t) по участкам монотонности. При таком выборе уравнение (3) будет представлять собой математическую модель реального временного ряда x,(t).

В работе детально проанализирован точно решаемый случай, когда (X (t) представляет собой Марковский дихотомический процесс (или D — шум), т.е. случайный процесс с двумя состояниями. В радиотехнической литературе его еще называют случайным телеграфным сигналом и определяют как случайную ступенчатую функцию, которая вслучайные моменты времени принимает значения cc(t) =±1 (или ±а , где а некоторая постоянная) с вероятностью 1/2. Скачки от одного значения к другому однородно распределены по времени и их число подчиняется статистике Пуассона, т.е. вероятность P(N) того, что за время dt произойдет ровно N скачков равна Р ( N ) = {у dt )'' е ""* / ¡V !. Параметр у характеризует среднюю частоту скачков в единицу времени и, с другой стороны, определяет характерное время т, ~ 'А спада корреляции процесса а (t) . Для рассматриваемого процесса его среднее <oc{i) >=0 и автокорреляционная функция < a(t + r)a(t) >= ехр(-и|/|), где операция о означает усреднение по ансамблю реализаций случайного процесса «(/).

Задача статистического описания динамической системы (3) при случайном а (I) заключается в определении вероятностных характеристик переменной x(t), характеризующей изменения состояний рассматриваемой экосистемы во времени. Для одноточечной плотности вероятности W (х, t) = ( w ( л, i ) ), где W (х, l)dx —вероятность того, что значение вектора состояния x(t) попадет в область (х, х + dx), а функция w{x,t) задает распределение переменной х в каждой отдельной реализации случайного процесса <2(7) и подчиняется стохастическому уравнению Лиувилля. Показано, что для искомого распределения W(x,t) система уравнений имеет вид:

\ д т . . 1 д ,+1 . ,

-+--+ +--(К'*' -К'-'Ш, = о

д1 2дх, ' ' 2 дх. ' '

-1 + уГР.+--(К( ) + Е" ш +--- к" ш = о

д1 1 2 . ' ' 1 2 '

/ /

где И7, =< а(1)\\>(х, /) > . Система (4) является точной и замкнутой. Полученная система кинетических уравнений является исходной для определения одноточечной плотности вероятностей состояний системы (3), представляющей из себя стохастическую смесь динамических процессов (1) и (2). Для рассматриваемой в работе модели О — шума могут быть найдены точные кинетические уравнения и для многоточечных распределений. Важным свойством реальных временных рядов является то, что они представляют собой одномерные функции времени. В этом случае система (4) существенно упрощается и описывает некоторый одномерный случайный процесс немарковского типа. Показано, что ее стационарное решение при произвольных нелинейностях Fc+>(x) и F(")(л;) имеет вид:

\FM-FH\ ( у W (x) = N . „ , exp --J

I 2

1 1

рм р(~)

\

dx

)

(5)

где N — нормировочная постоянная.

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

примера режима F(*' выбрана логистическая модель Ферхюльста •

х - ах - Ьхг, а режима F— экспоненциальное убывание численности х = -кх • Первое уравнение характеризует изменение численности при благоприятных условиях. Оно допускает экспоненциальный рост численности популяции на малых временах t и достижение некоторого стационарного значения хз = alb при достаточно больших значениях времени t. Второе уравнение описывает динамику численности популяции при неблагоприятных (катастрофических) условиях среды обитания, когда численность популяции экспоненциально убывает.

При перемешивании "благоприятных" и "неблагоприятных" условий по случайному закону Б — шума из (5) получается следующее стационарное распределение для численности популяции:

\а + к - Ьх\ \а - Ъх I

W|=N-!-!-7-, (6)

" 1« (*-«.) /2дА I , I1"- ' ^

Анализ (6) показывает, что статистическое перемешивание простых динамических режимов х - -кх и х = ах - Ьх2 приводит к большому разнообразию стационарных распределений для численности популяций. Возможные виды стационарных распределений при различных соотношениях между декрементом убывания численности (1/к) и инкрементом ее роста (1/а) и частотой стохастических переключений у рассмотренных динамик представлены на фазовой диаграмме (Рис. 3).

Рис. 3.

«Разовая плоскость и форма стационарных распределений численности популяции в

г<'> ; г(+) А 2

модели со стохастическими переключениями, '' = -кх , г = ах-ох ; правая

граница распределения определяется величиной х1 =а/Ь, /5 = к/с1, ц = у/2к-

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

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

В заключение главы обсуждается динамика изменения численности при периодическом чередовании более или менее благоприятных условий, соответствующих различным параметрам логистической модели. Если Fw(x,t) = ax-bx2, F{'\x,t) = px-qx1,a a(t) = ±1 является периодической функцией с периодом Q a(t) = a(t + 0) ),то система (3) имеет точные решения на полупериодах:

al

ае 1 Х(0 = *,(0 =-Г~1Г,если ' 6 [пв,(п + Цв]г

сх + Ье о >

pt

ре j

х(0 н х2(0 =-^ , если t е[(п+ -)9, (и + 1)0 ]

с2+дер 2

(7)

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

колебания численности популяции с периодом 0. Проведен подробный численный анализ задачи при различных соотношениях параметров. Показано, что как и в случае стохастического перемешивания, характер установившегося поведения (в данном случае амплитуда, частота и форма релаксационных колебаний, а также уровень, относительно которого происходят колебания) зависит от соотношения характерных времен Ир, На и частоты 2п ! в ■ Частота 2л I в является, таким образом, управляющим параметром и при периодическом перемешивании логистических динамик.

Рассмотрен частный случай, имеющий практическое значение, когда численность популяции в одном из режимов экспоненциально уменьшается со временем, а во втором ведет себя согласно логистической модели,т.е.если (х,*) = -кх, (х, 0 = ах-Ьх2.Показано, что существует предел, зависящий от характерных времен 1/к, 1/а, и частоты 2л I в , начиная с которого изменения численности могут стать необратимыми и привести к вымиранию популяции. На рисунке 4 приведен частный случай сценария изменения численности популяции при определенной частоте переключения динамик.

Кривые поведения численности популяции при периодическом чередовании благоприятных и неблагоприятных условий в зависимости от соотношения параметров а ! к .Горизонтальная асимптота х =а/Ь ,0 = 0.5 .1. а!к> 1.74 ; 2. я/*=1.74 ; 3. а / к <1.74 ;4. а/И<<1.74 .

Третья глава посвящена описанию разработанного автором программного комплекса "Система математического моделирования временных рядов" (СММВР). Комплекс предназначен для использования на машинах класса PC 386 и выше, как в среде DOS, так и WINDOWS.

В настоящее время пакет СММВР протестирован более чем на 50 временных рядах. При проведении расчетов получены достаточно устойчивые результаты. С помощью комплекса программ можно проводить моделирование временных рядов из различных областей знания. Результатом решения задачи являются аналитические модели, с некоторой точностью описывающие поведение данного временного ряда на выделенных промежутках. Алгоритмом программы предусмотрено построение экстраполяции временного рада на промежуток времени равный 1/4 выделенного при построении модели времени. На рисунке 5 приведен BP "Численность мохноногого курганника" с наложенной на него модельной кривой и участком экстраполяции. Участки роста модели описываются уравнением x(t) = с,ем", а участки убывания— x(t) = сге~ола>' Средний промежуток времени возрастания численности популяции 3.3 ± 1.1 года, средний промежуток убывания — 2.6 ± 0.55 года. Средняя относительная погрешность составляет 20.6%.

¡C/=20.V -.........Л------ ~ .... —. . .............. .............1-------------- " ..... ;

J

{

Г / !

S /

i ! ¡/г.

5 j i J

1 1 ^¡¡¡яШ _______1_______ ------j

) i T-.....1 . . j^oiSp^

Рис. 5.

Численность мохноногого курганника 1971 — 1996 гг. Тува (данные В.И. Забелина ТИКОПР СО РАН) Л — фактические данные, 2 — модель, 3 — участок экстраполяции.

Программы расчетов, входящие в состав комплекса, написаны на языке СИ 3.1, программа поддержки базы данных написана на языке CLIPPER 5.01. Пакет СММВР представлен множеством кодовых файлов: asmm00.exe; asmm01.exe; asmm02.exe; asmm03.exe; dbfedit.exe, листинги программ составляют более 50 страниц (около 2.5 тысяч

операторов). Программа asmm00.exe, выполняет функции диспетчера; программа asmm01.exe — поддерживает графические заставки пакета; программа asmm02.exe — реализует вычислительную часть пакета и построение модели, выбирая модельные кривые из следующих параметрических семейств: /(/,©) = ©0 +/в| — степенные функции и /(/,©) = ©0ев1' —экспоненциальные функции; asmm03.exe — программа поддержки системы помощи и сМе<Ш.ехе—редактор исходных файлов в формате "ёЬ?'.

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

Для обеспечения диалога с пользователем программа представляет следующие средства: систему меню для выбора альтернативных действий; систему запросов для выбора исходной информации и вида моделирующих функций; систему графического отображения области моделирования; систему помощи по эксплуатации программы. На рисунке 6 приведен рабочий экран программного комплекса СММВР.

(( Р Л П Л 'Л »1 VI г. г

тЩйМРгм(ММС* н

-¿и.

333 ;1Х<= О. ..... ...----•• -

;----- 1 1 -------- -------1-А2— -—.А ______г- А......:

Л; ------у? VI- —дД-Х- XV — /-(•---.¡Л \\ / /____ \---J/ щ Х-/.

Р \\/\ \ //

V ' * 1

313 ч/ ' 1 ' У "1 ...... и 1

е.а *-о.а|

4.6 «-0.2191 I ДМ

V». о.ичто х» <*» & а V*е -О 11С1 х1,за I (З4': { I * ™1

шЯШШШШШШШШЛ.

Рис б.

Рабочий экран программы. Отображены временной рад "Концентрация углекислого газа в атмосфере" и его модель. Т>% — относительная погрешность в процентах, ТТ, Т4 — средняя продолжительность участка роста, убывания соответственно, УГ, У4-—функции , ДГ=Л

Далее в главе приводится полученное аналитическое выражение для экстраполирующей функции:

хам) = х(о++ р*цм - +4^+,), (8)

где Х(/(+|) —значение ВР в /+/ точке, — значение ВР в /

точке, р+ — вероятность появления участка роста, Fí*\Fl') — аппроксимирующие функции, рассмотренные выше.

б

Рис 7.

а. Концентрация озона между 60° с.ш. и 69° ю.ш. по данным TOMS, б. Модель рада, приведенного на рис. а), наложенная на изображение рада.

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

чения BP и стоящих за ними процессов и явлений, т.е. представить временной ряд как результат перемешивания, чаще всего стохастического, двух простых динамик — возрастающей и убывающей. Положенный в основу моделирования BP анализ соотношений между характерными временами возрастания и убывания и частотой переключения режимов позволяет успешно выделять особенности поведения ряда, включая особенности вероятностных распределений. Например, на рисунке 7 изображен временной ряд "Концентрации озона" и его модель, отчетливо выделяется аномальное снижение содержания озона в 1982 году.

Развиваемый подход к моделированию BP позволяет достигнуть достаточно большой точности при относительной простоте построения модели, т. к. выбор моделирующих функций Fи/1'-1 возможениз класса простых гладких кривых, не имеющих особенностей, таких как степенные и экспоненциальные функции. Следует отметить, что рассмотрение временного ряда как некоторой суперпозиции двух динамик (возрастающей и убывающей), перемешиваемых по некоторому закону a{t), позволяет моделировать временные ряды любой степени сложности.

ВЫВОДЫ:

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

2. Представлены точно решаемые модели BP при стохастическом и периодическом перемешивании.

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

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

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

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

Положения, выносимые на защиту:

1. Для моделирования и экстраполяции временных рядов предлагается подход, основанный на представлении временного ряда как суперпозиции некоторых детерминированных и (или) стохастических монотонных состояний.

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

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

4. Программный комплекс, позволяющий на базе разрабаты-

ваемого подхода проводить математическое моделирование и экстраполяцию временных рядов.

ПО ТЕМЕ ДИССЕРТАЦИИ ОПУБЛИКОВАНЫ СЛЕДУЮЩИЕ РАБОТЫ:

1. Логинов В.М., Калуш Ю.А. Периодическое перемешивание динамических

поведений в экосистемах //Комплексное изучение аридной зоны Центральной Азии (Материалы международного рабочего совещания, Кызыл, 1994). Кызыл, 1998, с. 111-116.

2. Логинов В.М., Калуш Ю.А. Новый подход к математическому моделирова-

нию динамики экосистем //Сибирский экологический журнал, т.2, № 3,

1995, с. 196-210.

3. Логинов В.М., Калуш Ю.А. Математическое моделирование временных

рядов, возникающих при мониторинге природных процессов //Оптика атмосферы и океана, т.9, № 5,1996, с. 681-687.

4. Логинов В.М., Калуш Ю.А. Моделирование рядов долговременных

наблюдений //Труды 4 международного симпозиума «Эксперимент Убсу-Нур». М.: Интеллект, 1996. с. 127-135.

5. Логинов В.М., Калуш Ю.А. Программный комплекс для анализа и

прогнозирования временных рядов /ЛИ Межреспубликанский симпозиум

1996, с. 215-216.

6. Логинов В.М., Калуш Ю.А. Новый метод моделирования временных рядов

//Второй сибирский конгресс по прикладной и индустриальной математике (ИНПРИМ — 96), Новосибирск, 1996, с. 10-11.

7. Логинов В. М., Калуш Ю.А. Автоматизированная система математического

моделирования временных рядов природных процессов //Природные условия, история и культура Западной Монголии и сопредельных регионов (тезисы докладов III Международной научной конференции, Ховд, Монголия, 1997), Томск, 1997, с. 47-48.

8. А. с. 970532, Россия,/Система математического моделирования временных

рядов /В.М. Логинов, Ю.А. Калуш. 1997.

г

ТУВИНСКИЙ ИНСТИТУТ КОМПЛЕКСНОГО ОСВОЕНИЯ ПРИРОДНЫХ РЕСУРСОВ СО РАН

КАЛУШ ЮРИЙ АЛЕКСАНДРОВИЧ

НОВЫЙ ПОДХОД К МАТЕМАТИЧЕСКОМУ МОДЕЛИРОВАНИЮ ВРЕМЕННЫХ РЯДОВ В ЭКОЛОГИИ

05.13.16 —применение вычислительной техники, математического моделирования и математических методов в научных исследованиях (в экологии)

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

Кызыл— 1998

СОДЕРЖАНИЕ

1. Введение 3

1.1. Объект исследования 5

1.2. Краткая формулировка задачи и описание метода 6

1.3. Содержание работы 15

2. Обзор существующих методов исследования временных рядов 17

2.1. Анализ функции тренда 21

2.2. Анализ периодических (сезонных) колебаний 26

2.3. Анализ случайной составляющей 29

2.4. Обсуждение 30

3. Анализ и моделирование временных рядов 37

3.1. Математическая постановка задачи 3 7

3.2. Случайные переключения простых динамик 39

3.3. Случай периодических переключений простых динамик 52

3.4. Динамика изменения численности в условиях периодических чередований благоприятных и критически неблагоприятных условий 57

3.5. Обсуждение 60

4. Автоматизированная система моделирования временных рядов 64

4.1. Общее описание программного комплекса 64

4.2. Алгоритм построения модели 66

4.3. Структура программного комплекса «Система математического 69 моделирования временных рядов»

4.4. Функциональные характеристики программного комплекса СММВР 73

4.5. Аналитическое выражение для экстраполирующей функции 77

4.6. Обсуждение 81

5. Заключение 86

6. Литература 88

7. Приложение 1 98

1. ВВЕДЕНИЕ

Последние десятилетия свидетельствуют о возникновении своеобразного "экологического взрыва" — во всем мире резко возрос интерес к экологическим проблемам (см., например [1—15] и цитируемую там литературу). Приблизительно 20 лет назад экологии стали придавать значение, которое далеко выходит за рамки определения ее как раздела биологии. Изучая взаимосвязи между организмами и окружающей их средой, круговорот веществ и потоки энергии, делающих возможным существование жизни на Земле, экология связывает между собой различные области естествознания. Вполне естественно, что одни из самых мощных методов современного естествознания — математические методы, — стали широко применяться для решения экологических проблем. Возникла и бурно развивается наука — математическая экология [2, 4, 13, 16 — 19]. Одной из центральных проблем экологии вообще, и математической экологии в частности, является проблема устойчивости экосистем [20 — 24].

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

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

При математическом моделировании экосистем привлекается весь арсенал достижений современной математики, в частности, аппарат теории нелинейных дифференциальных уравнений в обыкновенных и частных производных, стохастическое исчисление, разнообразные теоретико-вероятностные подходы, численные методы и т.д. (см., например, [3,4,8—10,17,20,21,25—29]). Экспериментальной основой математических методов анализа экосистем являются наблюдения за поведением изучаемой экосистемы. Измеряемые характеристики экосистем, чаще всего протяженны по времени и носят как непрерывный, так и дискретный характер. Множество таких показателей, протяженных во времени называют временным рядом (ВР). Временные ряды являются важной частью изучения процессов, происходящих в экосистемах. Математические аспекты анализа ВР можно найти в [30 — 55].

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

1.1. Объект исследования

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

Временной ряд, полученный в результате наблюдения или эксперимента, представляет из себя некоторую кривую, отражающую изменение свойств процесса во времени. Для правильной интерпретации процесса, важно создание адекватной математической модели изучаемого ВР, к тому же, обладающей значительными предсказательными возможностями. Традиционный подход к моделированию ВР обычно заключается в следующем [32]:

• построение на основании имеющихся данных пар вида (х,у) некоторой кривой у = у(х), отображающей временной ряд,

• прослеживание характера полученной кривой, отслеживание экстремумов,

• изучение поведения условных средних значений на некоторых участках

группирования (выделение тренда).

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

1.2 Краткая формулировка задачи и описание метода

В данной работе формулируется новый универсальный подход к математическому моделированию временных рядов, понимаемых как совокупность наблюдений х (^ ), х (¿2),..., х(^„) исследуемой величины х(г), произведенных в последовательные моменты времени В основу этого

подхода положено несколько простых экспериментальных фактов, которые являются общими для всех ВР, возникающих при обработке данных мониторинга

различных природных процессов [56, 57]. Суть подхода проиллюстрируем на конкретных примерах временных рядов (Рис. 1.1).

j

1

1 V

щ

1_й 1 чЛ г. Л $ i

i A........ 1 Kt\ff J п . mг

1 ' г I \л { f-f м 1

П in ^ Р1 1 Ш \ тл / Ш* t T~ 4

П У\ У Л h /и/ г

f ' * * II_ " Ti

1995

М

г л

1.2

1963

¡д _ „

/ X \

" 7 / \ .. . L \ д

/ \ i ___/ .. т _i /

/ / i.. / \ X \s ч ....._______- /

к / Л*"""" \ / \

\ / X / \ / \ /

ч ..............V 1 \/

V

I960

Рис. 1.1

Примеры временных рядов. 1а) ВР «Глобальные изменения температуры» [8, 68]. 16) ВР «Концентрация углекислого газа в атмосфере» [8, 68]. 1в) ВР «Уровень реки Большой Енисей» (Республика Тыва). По данным гидрометеорологической

станции.

На рисунке 1.1а) приведено изменение глобальной температуры в последние несколько десятилетий. Не трудно заметить рост температуры, испытывающей случайные колебания. На рисунке 1.16) просматриваются детерминированные (периодические) колебания содержания углекислого газа в атмосфере, сопровождаемые общим ростом его содержания, а на рисунке 1.1 в) стохастические колебания уровня реки Большой Енисей в районе пункта наблюдения (с. Тора-Хем, Республика Тыва). Характерной особенностью приведенных временных рядов является то, что существуют промежутки времени, на которых изучаемая характеристика возрастает и существуют промежутки времени, на которых эта характеристика, наоборот, убывает. Закон смены участков возрастания и убывания может быть случайным, детерминированным, а также представлять их смесь. На участках монотонности временной ряд может быть описан некоторой непрерывной, гладкой кривой. Подход базируется на идее динамических и (или) стохастических переключений между режимами монотонности в точках смены знака производной [56, 57].

Суть предлагаемого метода обработки и моделирования BP удобно проиллюстрировать на основе схематического временного ряда, описывающего динамику поведения некоторого объекта или процесса {Рис. 1.2а) [57].

1.2а.

1.26.

1.2в.

Рис. 1.2.

а. Схематический временной ряд. б,в. "Веера лучей", отображающие периоды роста и убывания ВР.

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

рядов (рассмотренные ряды приведены в Приложении 1, для краткости П1).

— далее

Вероятно ста появления значения утаа раскрытия

I I ■!!■

I я

Угол р аскригия

Рис. 1.3.

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

Подсчет выполнен для рядов: «Глобальное изменение температуры» (П1, Рис. 1), «Концентрация озона в атмосфере» (П1, Рис. 2), «Содержание углекислого газа в атмосфере» (П1, Рис. 3), «Уровень реки Большой Енисей» (П1, Рис. 4), «Сейсмическая активность на территории Горного Алтая» (П1, Рис. 5), «Сейсмическая активность хребта Обручева» (П1, Рис. 6). Для реализации расчета выполнено преобразование подобия, заключающееся в том, что каждый ряд нормирован на свое максимальное значение, а минимальный участок роста (убывания) принят за единицу.

Из рисунков 1.26 и 1,2в видно, что для ВР всегда можно выбрать тем или иным способом некоторый эффективный (средний) угол наклона луча. Далее полагаем, что на участках возрастания ВР, изменение функции х(1) происходит по "закону эффективного луча". Обозначим получившуюся кусочно-гладкую

функцию через • Принимается, что эта функция описывает ВР х(1) на

10

участках возрастания. Повторяя рассуждения для промежутков времени, где ВР убывает {Рис. 1.2.в), получаем функцию /^(О, которая характеризует ВР х(1) с некоторой точностью на участках убывания. В результате ВР х(£) на временах наблюдения может быть апроксимирован некоторым эффективным временным рядом (0 + Р{ ) (/). Причем по построению, точность приближения зависит от того как выбираются "эффективные лучи". К достоинствам подобного построения можно отнести то, что для описания поведения ВР функции Р(+>(/) и могут быть представлены простыми гладкими кривыми, например, степенными (в частности линейными) и экспоненциальными функциями.

На рисунке 1.4 приведен пример временного ряда концентрации озона {Рис. 1.4а) и его модельное представление, наложенное на изображение ряда {Рис. 1.46). В качестве функции /7'ч,(7) выбрана функция /(/) = 12.2/096, а функции

— /{() = -7.94/105. Среднеквадратичная ошибка между реальным и модельным ВР составляет 1.2 % (изображение растянуто по вертикали из соображений наглядности). Это позволяет говорить о достаточной степени точности моделирования по предложенной схеме с использованием непрерывных гладких кривых. Как видно из рисунка, на модели четко разделяются три участка изменения концентрации озона: рост его содержания до 1982 года, резкое падение концентрации в 1982 году и постепенное замедляющееся восстановление содержания озона к прежнему уровню в последующие годы. Обсуждению вопросов, связанных с изменением концентрации озона в последние годы посвящен обзор [69].

озон

а.

1 и DИ= lTS

I А Д

. л КШ\1 \- - - - - -......- h Д/у \/!/ v \ ■ | У Vv ' \ ..«, s д л д д д Д,Д Д А Л А1 Л У\/V V V V-v-V V v-V v v V v О® 1 i 1 i 1 , t"

1976 1990

б.

Рис. 1.4.

а. Концентрация озона между 60° с.ш. и 69° ю.ш. по данным TOMS [69]. б. Модель ряда, приведенного на рис. а), наложенная на изображение ряда.

Примечательно, что не смотря на разнообразие проявлений и природу происхождения, характер поведения всех рассмотренных ВР может быть

охарактеризован поведением функций и а также законом

чередования режимов возрастания и убывания и их протяженностью. Введя индикаторную функцию а (/), принимающую значения ±1 в зависимости от

роста или убывания временного ряда, временной ряд можно уподобить некоторому дихотомическому процессу, характеризующемуся двумя состояниями: состоянием роста а = +1 и состоянием убывания а = -1. Участки монотонности временного ряда можно охарактеризовать некоторым средним показателем. Следует отметить, что статистика событий во времени, связанных с многими природными процессами с хорошим приближением может быть описана распределением Пуассона. Это распределение является предельным распределением для ряда дискретных распределений, таких как гипергеометрическое, биноминальное и так далее, когда число испытаний п велико. На Рис. 1.5 представлено распределение длин периодов временных рядов «Теоретическое распределение», ВР «Изменения глобальной температуры», ВР «Численность совы» (по данным [9]), ВР «Уровень реки Большой Енисей» в предположении, что среднее значение совпадает со средним участков монотонности.

5

о

X

н к о о.

О)

ш

0,8 0,6 0,4 0,2 0

НЕ

П-П

ИРяд1

■ Ряд2

□ РядЗ

□ Ряд4

2 3

Длина периода

Рис. 1.5

Распределения длин периодов при А.=1.

Ряд 1 - Теоретическое распределение. Ряд 2 —■ Распределение периодов ВР «Изменения глобальной температуры». Ряд 3 — Распределение периодов ВР «Численность совы». Ряд 4 — Уровень реки Большой Енисей. Как видно из приведенной диаграммы, распределения длин периодов реальных

ВР достаточно хорошо соответствуют теоретическому распределению, причем

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

Для экстраполяции ВР выбрана следующая схема: на экстраполируемом участке некоторым случайным образом строятся периоды возрастания и убывания ряда, подчиняющиеся распределению Пуассона со средними Г(+) (средний период возрастания) и Г(_) (средний период убывания). Поведение ряда на периодах роста и убывания описывается функциями Р<+> ([) и Р() (1). Полученная реализация ряда является случайной и не может индивидуально представлять его поведение в дальнейшем, поэтому строится некоторое множество реализаций по указанному сценарию, а затем производится усреднение по ансамблю реализаций. Указанное усреднение представляет собой экстраполяцию поведения временного ряда. Расчеты показали, что ошибка экстраполяции со