Применение КА-модели для исследования влияния загрязнений на динамику популяций голомянок и макрогектопуса в озере Байкал Текст научной статьи по специальности «Математика»

2014 Дискретные модели реальных процессов №1(23)

ДИСКРЕТНЫЕ МОДЕЛИ РЕАЛЬНЫХ ПРОЦЕССОВ

УДК 519.17

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

И. В. Афанасьев

Институт вычислительной математики и математической геофизики, г. Новосибирск,

Россия

E-mail: ivafanas@gmail.com

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

Ключевые слова: самоорганизация, дискретное моделирование, клеточный автомат, композиционные модели, хищник — жертва.

Введение

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

В [4] предложена клеточно-автоматная модель (КА-модель) динамики численности восьми групп организмов озера Байкал, позволяющая учитывать пространственное распределение особей и возможные локальные загрязнения. В этой модели область моделирования является квадратной, и не учитываются влияние водных течений и сезонные изменения.

В данной работе КА-модель [4] модифицирована. Добавлено соответствие области моделирования пространственной структуре озера Байкал, учтено влияние сезонных изменений на динамику численности и перемещений организмов под действием водных течений. Скорректирован метод учёта влияния локальных загрязнений. Представлены результаты верификации модели. Приведены оценки минимальной интенсивности

1 Исследование выполнено по Программе фундаментальных исследований Президиума РАН, проект № 15.9 (2012г.), а также при частичной финансовой поддержке РФФИ в рамках научного проекта №11-01-00567а.

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

1. Композиционная КА-модель

Рассмотрим три вида организмов: макрогектопус (шасгоЬе^ориз), малую (сошерЬо-гия dybovski) и большую (сошерЬогиз Ьа1ка1еп81з) голомянок.

Каждый из видов разделён на возрастные группы (вид будем обозначать буквой из |т,^,6}, возрастную группу — цифрой из {1, 2, 3}):

— макрогектопус — неполовозрелые Ш1, половозрелые т2;

— малая голомянка — однолетки ^1, неполовозрелые ^2, половозрелые ^3;

— большая голомянка — однолетки Ь1, неполовозрелые Ь2, половозрелые Ь3.

Между группами определены взаимоотношения хищник — жертва и демографические взаимоотношения.

КА-модель динамики популяций организмов озера Байкал определяется четвёркой понятий

N = (Е,М,^,р),

где Е — алфавит состояний клеток; М — множество имён клеток; ^ — глобальный оператор перехода; р — режим функционирования.

Модель — параллельная композиция [5] восьми КА, каждый из которых предназначен для моделирования численности конкретной группы организмов.

Пусть Q — квадратная сетка, покрывающая озеро Байкал. Множество имён клеток М — объединение восьми попарно непересекающихся множеств М“:

М = МГ и М2Г и М^ и М2 и м3 и МЬ и м2 и М3Ь.

Предполагается, что каждой ячейке из Q соответствует восемь имён клеток — по одной из каждого из множеств М“, т. е. для всех г и а существует биекция ф“ между множеством М" и множеством ячеек сетки Q:

фа : Q ^ ма.

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

5(с) = ((^1(0), П1),..., (с),щ))

называется локальной конфигурацией, где щ Е Е — состояния клеток; ^ : М ^ М — именующие функции, т. е. функции, определяющие для клетки с именем с имена клеток, с ней взаимодействующих.

В общем случае локальный оператор перехода / : {5(с)} ^ {5(с)} принимает на вход локальную конфигурацию Б (с). Его результат — локальная конфигурация / (Б (с)), в которой изменены состояния клеток (возможно, не всех) из Б (с). Итерация, или применение глобального оператора ^ — применение локального оператора / ко всем клеткам. В модели глобальный оператор ^ является последовательной композицией двух операторов [5]

^ о ^2,

где предназначен для моделирования перемещения организмов, ^2 —для моделирования процессов поедания, вымирания и роста. Оператор перемещения ^\, в свою

очередь, является последовательной композицией оператора собственного перемещения Fd и оператора перемещения под действием течений /,:

^ о /,.

Известны два основных режима применения глобального оператора: синхронный и асинхронный. Синхронный режим предполагает, что сначала вычисляются новые состояния клеток согласно локальному оператору перехода, а затем все клетки одновременно изменяют свои старые состояния на новые. В асинхронном режиме случайно выбирается клетка, вычисляется ее новое состояние и старое состояние клетки сразу заменяется на новое [6].

В предложенной КА-модели используются оба режима. Сначала асинхронно применяется /\, затем синхронно применяется /2.

1.1. Собственное перемещение организмов

Пусть /г и ^ —локальный и глобальный операторы целочисленной диффузии, предложенной в работе [7],

/г : {Б1(с)} ^ {Б1(с)}

где Б1(с) —набор соседей клетки с именем с, включая саму эту клетку. Соседями назовем клетки с именами с“,с“/ Е М“, соответствующими соседним ячейкам сетки Q. Клетки с именами из разных подмножеств М“ соседями не являются, т. е. если с“ Е М“ и св Е Мв — соседи, то г = ], а = в.

Применение /г к клетке (с, п) выполняется по следующему алгоритму:

1) Пусть (с1, щ1),... , (ск, пк) — соседи клетки (с, п), к ^ 4.

Случайно равновероятно выбирается (с,щ).

2) Новые состояния п/ и п^ вычисляются по формулам

п/ = п — [ап] + [ап ], Щ = щ + [ап] — [ап ],

где а — коэффициент целочисленной диффузии, 0 ^ а ^ 1.

Оператор ^ применяется в асинхронном режиме.

Оператор определяется следующим образом. Пусть I — физический размер квадратной ячейки, гсг — крейсерская скорость организмов вида а возраста г, Д£ — физическое время, соответствующее одной глобальной итерации КА. Тогда максимальное число ячеек, пройденное организмом за время Д£, может быть вычислено как

гсгД

К‘ = “Г'

Так как применение /г перемещает организм только на одну клетку, то перемещение с крейсерской скоростью достигается К“-кратным применением оператора на множестве М“:

^|м« = (^)К“ '

1.2. Перемещение под действием течений

Пусть з£геат(с) : М ^ К2 — карта течений — отображение множества имён клеток М в пространство двумерных вещественных векторов.

Вектор з£геат(с) = (гх,гу) — физическая скорость течения в клетке с именем с; х соответствует направлению запад — восток, направление у — юг — север.

Локальный оператор /, моделирует перемещение организмов вдоль траектории течения. Алгоритм применения /, к клетке с именем с рекурсивен. Сначала просматривается клетка с именем с, затем её соседи, в которые течение переносит организмы, затем соседи соседей и так далее. Рекурсия обрывается на расстоянии, дальше которого течение не может перенести организмы за время, соответствующее одной итерации.

Пусть I — физический размер ячейки, Д* — физическое время, соответствующее одной итерации.

Алгоритм применения /, к клетке с именем с:

1) Найти сХ и ск — имена соседей клетки с именем ск в направлении гХ и , где к — уровень вложенности рекурсии; ск — имя просматриваемой клетки. На первом уровне рекурсии просматривается клетка с именем с (с1 = с). На втором уровне просматриваются две клетки — соседи клетки с именем с1 , в которые попадает течение (с2 Е {сХ, суУ}). На к-м уровне просматривается не более 2к клеток, (гХ, гк) = з£геат(ск) —вектор течения в клетке с именем ск.

2) Вычислить пХ, пк — число организмов клетки с именем ск, перемещаемых тече-

к к

нием в клетки с именами сХ и с^ соответственно:

гк

Iпк | к | ХI к I, если «Х >1,

I 1 гХ1 + 1 гк 1 п

Л гк пУ

Х

к

п\ к, Г, к,^ если«Х ^ 1, 1 гХ1 + 1 гк 1

гк

к У к 1

п\------:-------:---г~Г, если «к > 1,

I гХ | + | гк Г у ’

гк

к У к к ^ 1

п^П Г-:—^ Г-Г «к, если ^ 1,

1 гХ 1 + 1 гк 1

у,

У '

где пк — модельное число перемещаемых организмов клетки с именем ск. Число пк при к > 1 вычисляется на (к — 1)-м уровне рекурсии на шаге 3; п1 — состояние клетки с именем с, т. е. на первом уровне вложенности под действием течения перемещаются все особи, находящиеся в клетке с именем с; ^Х, «к —

к

модельные расстояния, преодолеваемые течением за время тк, в количестве клеток:

ткгк ткгк

„к = т гХ „к = ' гу.

«х I , „у I .

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

т1 = Д*, тк = тк-1 — Дтк-1;

Дтк-1 — физическое время, затраченное течением на преодоление клетки на (к — 1)-м уровне вложенности рекурсии:

Дтк-1 = ——- или Дтк-1 —

гХ-1 гк-1'

3) Если ^Х > 1 (т. е. течение за время тк может пересечь более одной клетки в направлении х), то алгоритм рекурсивно выполняется для клетки ск+1 = сХ с числом перемещаемых организмов пк+1 = пХ и временем, затраченным на перемещение, Дтк = —т. Аналогичная проверка выполняется в направлении у.

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

1.3. Оператор изменения численности Пусть /2 — локальный оператор изменения численности:

/2 : {^2(с)} ^ {Б2(с)},

где $2(с) —набор клеток-близнецов (рис. 1). Клетки с именами с? Е М? и св Е Мв

з

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

3 )-1(св:

(ф“) 1(с“) = (фв) 1(св). Таким образом, для с Е М? и 9 = (ф“) 1(с) Е ф

^2(с) = (ФГ(9),^т(9),^(9),^(9),^(9),^Ь (9),^Ь(9),^Ь (9)).

Рис. 1. Клетки-близнецы 52(с)

Локальный оператор /2 применяется синхронно. Новое состояние клетки с именем с" Е М? вычисляется по формуле

(п?)' = п“ + (р?п? — А“п“ — б? п?)Д*,

где ^ —возрастная группа особей, порождающих особей возраста г; Д* — физическое время, соответствующее одной итерации КА; р?па —приток в группу за счёт рождаемости или старения предыдущей группы; А"п" —отток из группы за счёт смертности; ^п? — отток из группы за счёт старения.

Коэффициенты рождаемости р?, смертности А? и старения взяты из работы [3]:

1) коэффициенты старения 0? постоянны;

2) коэффициенты смертности А2, А^, А2, А3 постоянны;

3) коэффициенты смертности Ат, Ат, А^, А? имеют вид

А? = а? + 6?(п2 + п3) + ^?(п2 + п3), (1)

где а?, 6", # —постоянные величины; 6? и характеризуют смертность от выедания хищниками (малыми и большими голомянками) и зависят от вкусовых предпочтений хищников; а? характеризует смертность от остальных причин;

4) коэффициент рождаемости рт постоянен;

5) коэффициенты рождаемости р1,р1 имеют вид

р? = ^>т+пт) + V “(п!+п1),

где постоянные коэффициенты и V? зависят от рациона хищников.

Чтобы учесть зависимость рождаемости от сезонов, в КА-модели коэффициенты рождаемости р^3 и р1 умножаются на периодические функции зеазопз : К ^ К и зеазопь : К ^ К, графики которых приведены на рис. 2. Период функций 8вазоп3(*) и звазопь(*) равен одному году: звазопь(* + 1) = звазопь(*), звазоп3(* + 1) = звазоп3(*).

2.5

2

1.5 1

0,5 0

О СО ЦП со со со ио со ю со цп со о) со

О т- СЧ 0~ СО хг Ш_ 0' Ю г- СО 0" О)

о о* о о* о’ о" о о" о' о"

Период, год

Рис. 2. Графики функций ввавощ(£) (кр. 1) и ввавоп,(^ (кр. 2). По оси абсцисс отложен один год. Ноль соответствует 1 января

2. Верификация модели

Верификация проводится по трём видам организмов (суммарно по возрастным группам) и следующим критериям:

— —а/Ва — отношение годовой продукции к среднегодовой биомассе организма а. Годовая продукция Ра — положительный прирост биомассы организма вида а £ {т, Ь} за год. В это значение входят массы родившихся организмов и прирост массы организмом во время взросления;

— Nа/Ыр — отношение числа организмов вида а к числу организмов вида в; Nа — усредненная по пространству среднегодовая численность;

рт=—т+—т > Вт = вт+вт , ыт=,

—, = Р1 + Р2 + рз, В, = В1 + В\ + В3, = N + N2 + N3,

Рь = —1 + —2 + —3, Вь = В-1 + В-2 + В-3, N = N1 + N2 + N3.

В [8] дана оценка биомассы макрогектопуса Вт в 110000 тонн. Оценка продукции макрогектопуса Рт варьируется от 330000 до 900000 тонн [8], т. е. коэффициент Рт/Вт должен находиться в пределах 3-8.

Оценка коэффициента ——- взята из [9] и равна 1,24. Оценки коэффициентов

N N В, + В-

—т и —т взяты из [3]. Результаты верификации представлены в таблице. Модельные

Nd Щ

оценки отличаются от оценок, приведённых в [3, 8, 9], не более чем на 20%.

Результаты верификации модели

Оценки Pm вт Pd + Pb Bd + Вь Nm Nd Nm Nb

В [3, 8, 9] 3-8 1,24 6,05 21,52

Модель 5,77 1,49 6,00 20,45

3. Влияние загрязнений

Модель позволяет исследовать влияние вероятного загрязнения.

Пусть poll(c) : M ^ R+ — карта загрязнений — функционал, ставящий в соответствие имени клетки положительное число, характеризующее интенсивность загрязнения. Предполагается, что загрязнение влияет на смертность организмов.

Новые коэффициенты смертности для хищников Ad, Adf,A2, A3 вычисляются по формуле

(A“ (c))' = A“ + A?poll(c).

Новые коэффициенты смертности для жертв Лт,Ат,А1,Л1 вычисляются по формуле

(Ла(с))/ = Ла + aapoll(c),

где а“ — постоянный коэффициент смертности из (1).

Разделение для коэффициентов на хищников Л^,Л^,Л2, Л3 и жертв А^А^А^АЗ при учёте загрязнений сделано для того, чтобы в множитель смертности от загрязнения не попала смертность от выедания хищниками из (1).

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

Рис. 3. Карта загрязнений: серый цвет — земля; белый — вода; более тёмный цвет означает большую интенсивность загрязнения

4. Применение модели для исследования влияния загрязнений

Для исследования влияния локальных загрязнений проведён вычислительный эксперимент. Начальное состояние — равномерное распределение особей по области моделирования, совпадающее по значению с устойчивым состоянием численностей, взятому из [3]. Значение ро11(с0) в эксперименте равно 14.

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

У Я Я

t= 50 t = 200 Г =2000

Рис. 4. Несколько итераций для половозрелого макрогектопуса. Более тёмный цвет означает большую плотность организмов

Динамика популяций в клетке на севере озера (рис. 5) стремится к устойчивому колебательному процессу с периодом колебаний 1 год. Годовые колебания — следствие сезонной зависимости коэффициентов рождаемости больших и малых голомянок.

£ ■

Ed -

5 -

с -

я : ]: -

I ::: _ § -

О

^ ::: -

Рис. 5. Динамика модельной плотности в клетке на севере озера (без загрязнения) для неполовозрелого макрогектопуса (а) и неполовозрелой малой голомянки (б)

Динамика популяций в эпицентре загрязнения представлена на рис. 6. Значение ро11(с0) загрязнения в эпицентре равно 14. Этого достаточно для того, чтобы все организмы погибли.

СМ ом со со со

Гэды

Рис. 6. Динамика модельной плотности в клетке в эпицентре загрязнения (кр. 1) для неполовозрелого макрогектопуса (а) и неполовозрелой малой голомянки (б); кр. 2 — аналогичная динамика в клетке без загрязнения

Динамика популяций в области загрязнения в клетке со значением интенсивности poll = 2,3 представлена на рис. 7. Численность жертв увеличилась, а численность хищников уменьшилась по сравнению со среднегодовым значением в области, где загрязнение не оказывает влияния.

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

Распределение интенсивности загрязнения и изменение численности макрогектопуса на 2000-й итерации алгоритма вдоль отрезка с юго-востока на северо-запад представлены на рис. 9.

По результатам экспериментов можно сделать следующие выводы:

— при poll (с) > 10 плотность макрогектопуса на 2000-й итерации падает ниже 3% от средней плотности в незагрязненной области, что позволяет говорить о практически полном вымирании организмов в загрязненной области;

— если poll (с) £ (0,15, 5), то в этой области увеличено число жертв и уменьшено число хищников; степень различий зависит от значения poll(m);

НО CD CD

CNI CM W fO

Годы

Годы

СМ СМ СО СО

Годы

5

ч

о

Годы

б

а

Рис. 7. Динамика модельной плотности в клетке со значением загрязнения poll = 2,3 (кр. 1) для неполовозрелого макрогектопуса (a) и неполовозрелой малой голомянки (б); кр. 2 — аналогичная динамика в клетке без загрязнения

Рис. 8. Отрезок, вдоль которого исследовалась динамика модельной плотности

— если ро11(с) < 0,03, то влияние загрязнения незаметно на фоне среднегодовых колебаний численности организмов.

Расстояние вдоль отрезка, км

Рис. 9. Распределение интенсивности загрязнения вдоль исследуемого отрезка (а) и модельной плотности неполовозрелого макрогетопуса на 2000-й итерации (б). По оси абсцисс отложено расстояние на отрезке от эпицентра загрязнения, правая граница соответствует 11,5 км

Заключение

Предложенная ранее КА-модель динамики численности восьми групп организмов озера Байкал скорректирована и дополнена учётом пространственных особенностей области моделирования, зависимости взаимодействий между группами от времени года и перемещения особей под действием течений.

Модель верифицирована. Результаты верификации по критериям отношения продукции к биомассе и частот встречаемости отличаются от оценок, приведённых в работах [3, 8, 9], не более чем на 20%.

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

ЛИТЕРАТУРА

1. Свирежев Ю. М., Логофет Д. О. Устойчивость биологических сообществ. М.: Наука, 1978. 352 с.

2. Базыкин А. Д. Нелинейная динамика взаимодействующих популяций. Ижевск: ИКИ, 2003. 368 с.

3. Зоркальцев В. И., Казазаева А. В., Мокрый И. В Модель взаимодействия трех пелагических видов организмов озера Байкал // Современные технологии. Системный анализ. Моделирование. 2008 №1. С. 182-193.

4. Афанасьев И. В. Клеточно-автоматная модель динамики численности организмов озера Байкал // Прикладная дискретная математика. 2012. №1. С. 121-132.

5. Bandman O. L. Cellular automata composition techniques for spatial dynamics simulation // Simulating Complex Systems by Cellular Automata. Understanding complex Systems / eds. A. G. Hoekstra et al. Berlin, 2010. P. 81-115.

6. Бандман О. Л. Клеточно-автоматные модели пространственной динамики // Системная информатика. 2006. №10. С. 58-113.

7. Medvedev Y. G. Multi-particle cellular automata models for diffusion simulation // Meth. Tools Parall. Program. Multicomp. 2011. V. 6083/2011. P. 204-211.

8. Мазепова Г. Ф, Тимошкин О. А., Мельник Н. Г. и др. Атлас и определитель пелагобионтов Байкала. Новосибирск: Наука, 1995. 693c.

9. Стариков Г. В. Голомянки Байкала. Новосибирск: Наука, 1977. 94 с.