Многокомпонентная модель для экосистемы Охотского моря Текст научной статьи по специальности «Математика»

Многокомпонентная модель экосистемы Охотского моря

Ишмукова И,В, ([email protected])

Тихоокеанский научно-исследовательский рыбохозяйственный центр (ТИНРО-центр)

Введение

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

Процесс построения и использования математических моделей включает в себя несколько этапов [19]. Хотя к настоящему моменту известно некоторое количество формальных моделей [1, 5, 6, 7, 9, 10, 16], специалисты по морской биологии в основном выражают свои представления об экосистемах в неформальном виде [13], [14], [15], К примерам такого рода относится и модель из [17], описывающая макроэкосистему Охотского моря, В этой модели кроме общей блок-схемы приведены значения выеденной биомассы соответствующими трофическими блоками в период 1980-х годов, Эти значения представляют собой "средний" установившийся поток масс, В математических терминах это означает, что при приведенных значениях в экосистеме установлено равновесное состояние и малые изменения фазовых переменных приводят к затухающим динамическим процессам и возврату в положение равновесия, Таким образом, проверка "корректности" этой неформальной модели должна включать в себя построение математической модели и исследование устойчивости точек равновесия.

Построение математической модели по известным из [3] и [17] данным включает в себя:

• выбор типов математических соотношений, описывающих обменные процессы;

В работе для формального описания обменных процессов используются известные способы описания динамики биомасс видов в трофических цепях [12], В каждом дифференциальном уравнении, соответствующем трофическому блоку, в правой его части описано потребление этим блоком органического вещества (составляющие со знаком " + "), выедание биомассы этого блока другими блоками, естественная и промысловая смертность (составляющие со знаком " —"), Трофические функции и функции смертности - линейные. Параметрическая идентификация проведена с точностью до 2-х параметров. Исследована локальная устойчивость точек равновесия. Рассмотрены динамические процессы, возникающие в системе при малом возмущении стационарных точек.

1 Описание моделей

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

Рис, 1, Потоки вещества в модели экосистемы Охотского моря в 80-е гг.

Пусть х0 - биомасса фитопланктона, х1 - биомасса бактериопланктона, х2 - биомасса простейших, хз - биомасса нехищного зоопланктона, х4 - биомасса хищного зоопланктона, х5 - биомасса пек топа, х6 - биомасса нехищного зообентоса, х7 - биомасса хищного зообентоса, х8 - биомасса нектобентоса, х9 - биомасса детрита, хю -биомасса биогенов. Предположим, что биосистема моделируется автономной системой обыкновенных дифференциальных уравнений, описывающей перенос вещества между блоками, вида

х = / (х), (1)

где х = (х0,..., хю) - биомассы видов в сообщеетве, / = (/0,..., /ю).

Рассмотрим две модели. Первая описывает тот случай, когда нет промысла, т.е. биомасса видов уменьшается только по причине естественной смертности, и правая часть /(х) = /(1)(х), где

/¿1)(х) = Р0(хю)х0 - ^ К0.(ж0)ж,- - т0(ж0)ж0,

/(1)(х) = ^ к.дУ^ж.)Х1 + ¿4,9^0(хд)ж1 - ^ У^(ж - Ш^ж^ж,

/¿^(ж) = ^ .У'г (ж.)Хг - ^ У,.(ж^ж. - Ш^ж^, 2 = 2...8,

8 8

/91)(Х) = ^ ^(1 - - ^ Уд. ^ж. - ^0(^9)^1 + ^ Ш^ж^,

г=1 ¿е^ г=0

/1(01)(х) = (1 - к1,9)^0(Х9)Х1 - Р0(хю)х0.

Здесь функция Р0(ж10) описывает скорость потребления биогенов (ж10) фитопланктоном (ж0) при фотосинтезе; фупкция ^0(ж9) описывает скорость потребления детрита (ж9) бактериоплапктоном (ж^; функция У. (жг) - функция выедания, описывающая скорость потребления биомассы г -го вида единицей биомассы j-тo вида; функция шг(жг) - скорость отмирания вида по причинам, не связанным с трофическими отношениями; - коэффициенты использования усвоенной пищи па рост, - множество номеров видов, которые выедаются г-ым видом; - множество номеров видов, которые выедают г-ый вид, г = 0,10 ] € 3г,

Во второй модели будем учитывать вылов основных видов. Важные для промысла объекты относятся к таким трофическим блокам, как нектон (ж5) и нектобентое (ж8). Пусть функции /г(2)(ж) = /г(1)(ж) для г = 0,4, 6, 7, 9,10, а

/б2)(х) = ^ к.,5У,-5(ж.)ж5 - ^ К5,.(ж5- тб(жб)жб - и5(ж5)ж5,

/8(ж) = ^ к.,8У.-,8(ж.)Ж8 - ^ Уу(Ж8- Ш8(ж8)ж8 - Цв^)^.

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

Р0 (ж 10) = Р0Ж10, ^0(^9) = 50^9,

шг (жг) = шг,

ига(жга) = ura,

где р0, д0, ггш^ и - положительные константы, г = 0,..., 9, € 3", I = 0,..., 8, п = 5, 8, Таким образом, сформированы 2 модели переноса вещества между блоками (см. (1) и (2))

ж = /(к)(ж), к = 1, 2. (3)

2 Идентификация моделей

После линеаризации функций из первой и второй модели, возникает задача параметрической идентификации. Проведем ее на основе данных, известных из [3] и [17].

Если х*(к) = (х1(к) , ...,хП(к)), к = 1, 2 характеризует состояние равновесия, то выполняется условие

/(к)(х*) = 0, к =1, 2. (4)

Из [3] и [17] известны средние "устоявшиеся" биомассы. Предположим, что они

соответствуют равновесным состояниям х*(к) = (х1(к),..., хП(к)), к = 1, 2, Коэффици-

(к)

енты т\- определяются из следующих предположении:

гЙЧ*(к)ж*(к) = р(к), * = 0,..., 9, 3 е У, к =1, 2, (5)

где р(к) - это величины выедания биомассы одного трофического блока другим (метки дуг графа на соответствующих рисунках). Определив коэффициенты т( вычисляем, используя условие (4) для /(1), коэффициенты ш(1), х^ и х^ с точностью до параметров р0 и I = 0,..., 8,

Для второй модели введем дополнительные уравнения для определения коэффициентов м5 и м8 :

ЩХп = Wn, (6)

где Wn - известные величины выловов (значения которых указаны на рис, 1), п = 5,8, Определив коэффициенты т( 2), м5 и м8, вычисляем, используя (4) для /(2), коэффициенты естественной смертности т(2), а также биомассу детрита х92) и биомассу биогенов ж 10 с точностью до р0 и д0, * = 0,..., 9, 3 е У , I = 0,..., 8,

Проведенные таким путем рассуждения позволяют вычислить все параметры моделей за исключением р0 и д0.

3 Исследование моделей

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

Предположим, что биомассы из [3] и [17] соответствуют стабильным значениям ж(к), к = 1, 2 в течение периода времени [0,Т], Исследуем устойчивость точек равновесия ж(к) методом линейных приближений [11], Пусть А(к)(ж(к)) - это спектр

матрицы А(к)(ж(к)) = /Хк)(ж), т.е. множество ее собственных значений. Если

УА(к) е £р А(к)(ж(к)) Де А(к) < 0, (7)

ж(к)

ЗА(к) е £р А(к)(ж(к)) Де А(к) > 0, к =1,2, (8)

то равновесие неустойчиво,

В табл. 1 приведены действительные части собственных чисел матрицы А(к)(ж(к)), к 1, 2, которые были получены при р0 и д0, равными 0.2 и 0.03 соответственно. Некоторые

Таблица 1, Действительные части собственных чисел

Де А(1) Де А(2)

-1.9180 -1.9622

-1.9180 -1.9622

-4.5521 -4.2834

-4.5521 -4.2834

0.3781 0.3786

0.3781 0.3786

-0.0373 -0.0354

-0.0373 -0.0354

0.0311 0.0287

0.0311 0.0287

1.6564 * 10-15 0.0009

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

ж(1) ж(2)

Однако даже если положение равновесия не является устойчивым, имеет смысл рассмотреть возникающие в системе динамические процессы. Интерес представляют решения, обладающие определенными свойствами (периодические, асимптотические, показатели Ляпунова, аттракторы и фракталы). Для этого численными методами были найдены решения задачи Коши при начальных условиях ж(0) = (1 + $)ж*(к), к = 1, 2. Рассмотрены различные варианты решений при изменении р0 и 90 от нуля до единицы,

В табл. 2 и в табл. 3 собраны значения параметров, когда р = 0.2, а 9 = 0.03,

где В - это биомасса (млн.т), К2 - коэффициент использования пищи па рост, Р

т(к)

коэффициенты выедания одного блока другим, т(к) - коэффициенты смертности, и - коэффициент уловиетоети, к = 1, 2,

При численном моделировании использовались нормированные данные. Нормировка проводилась следующим образом. Пусть ж*(к) - точка равновесия, к = 1, 2, Введем замену переменных

ж

^ =

ж

*(к) '

Таблица 2, Значения параметров для первой модели (без учета промысла)

Объекты исследований т(1) В К2 Р г(1) т(1)

Фитопланктон т(1) х0 61.860 - - - 54.61

Бактерии т(1) 1 28.58 к9)1 = 0.32 р9>1 = 20409 - 73.48

Простейшие т(1) т2 11.431 ко,2 = 0.55 Ро,2 = 4979 Го,2 = 7.04 242.17

к1>2 = 0.55 р1>2 = 2325 Г1,2 = 7.12

Нехищный зооплантон т(1) т3 314 ко,з = 0.4 Ро,з = 5521 Го,з = 0.28 8.35

к1,з = 0.4 Р1,з = 1932 Г1,з = 0.22

к2,з = 0.4 Р2,з = 1073 Г2,з = 0.299

кд,з = 0.4 Р9,з = 1379 гд,з = 0.0006

Хищный зоопланктон т(1) т4 115 км = 0.35 Р1,4 = 15 гм = 0.0046 3.007

кз,4 = 0.35 Рз,4 = 1050 гз,4 = 0.029

Нектон т(1) т5 33.5 кз,5 = 0.3 Рз,5 = 61 гз,5 = 0.006 0.69

к4,5 = 0.3 Р4,5 = 26 Г4,5 = 0.007

Нехищный 30- обен-тос т(1) т6 208.6 к2,б = 0.35 Р1,6 = 159 Г1,6 = 0.027 0.48

к2,б = 0.35 Р2,6 = 176 Г2,6 = 0.074

кз,б = 0.35 рз,6 = 227 гз,6 = 0.003

к9;6 = 0.35 Р9,6 = 196 г9;6 = 0.00013

Хищный зообентос т(1) 21.4 кб,7 = 0.3 Р6,7 = 147 Г6,7 = 0.033 2.23

Нектобентос т(1) т8 5 кз,8 = 0.2 Рз,8 = 1 гз>8 = 0.0006 1.072

к4,8 = 0.2 Р4,8 = 1 Г4,8 = 0.0017

к5,8 = 0.2 Р5,8 = 2.9 Г5,8 = 0.0173

кб,8 = 0.2 Р6,8 = 18.2 Г6,8 = 0.0174

к7,8 = 0.2 Р7,8 = 3.7 г7>8 = 0.0346

Детрит т(1) т9 7141.26 - - - -

Биогены т(1) т10 280.43 - - - -

Тогда системы (3) перепишутся в следующем виде

ги = Т^И, к = 1, 2. (9)

На рис, 2 изображены изменения переменных и>з, и>5 и и>8 во времени, а также фазовые портреты в плоскостях этих переменных. Функция и>3 при отклонении от точки равновесия 6 = 0.1 возрастает с небольшими колебаниями. Переменные и>5 и и>8 монотонно убывают при 6 = -0.1,

На рис, 3 изображены динамика переменных и>3, и>5 и и>8. Там же представлены фазовые портреты этих переменных. Как видно, при возмущении 6 = 0.1, функции монотонно возрастают, а при 6 = -0.1 - монотонно убывают. На графиках видно, что точки равновесия в первой и второй моделях неустойчивы.

Графики носят иллюстративный характер. Можно приводить достаточно большое количество вариаций подобных расчетов.

Таблица 3, Значения параметров для второй модели (е учетом промысла)

Объекты исследований x(2) В K2 P r(2) m(2) u

Бактерии x(2) Jb i 28.58 fc9>1 = 0.32 p9,1 = 20409 - 73.48 -

Простейшие x(2) x2 11.431 ko,2 = 0.55 Po,2 = 4979 ro,2 = 7.04 242.17 -

k1>2 = 0.55 p1,2 = 2325 Г1,2 = 7.12

Нехищный зоопланктон x(2) x3 314 ko,3 = 0.4 Po,3 = 5521 ro,3 = 0.284 8.354 -

ki,3 = 0.4 P1,3 = 1932 r1,3 = 0.22

k2,3 = 0.4 P2,3 = 1073 r2,3 = 0.299

k9,3 = 0.4 P9,3 = 1379 r9,3 = 0.0006

Хищный зоопланктон x(2) x4 115 ki,4 = 0.35 P1,4 = 15 r1,4 = 0.005 2.543 -

k3,4 = 0.35 P3,4 = 1050 r3,4 = 0.029

Нектон x(2) x5 33.5 &3,5 = 0.3 P3,5 = 61 r3,5 = 0.006 0.637 0.06

k4,5 = 0.3 P4,5 = 26 r4,5 = 0.007

Нехищный зообентос x(2) x6 208.6 ki,6 = 0.35 P1,6 = 159 r1,6 = 0.027 0.48

fc2,6 = 0.35 P2,6 =176 r2,6 = 0.074

^3,6 = 0.35 P3,6 = 227 r3,6 = 0.003

k9,6 = 0.35 P9,6 = 196 r9,6 = 0.00013

Хищный зообентос x(2) Jb 7 21.4 k6,7 = 0.35 P6,7 = 147 r6,7 = 0.033 2.23

Нектобентос x(2) x8 5 k3,8 = 0.2 P3,8 = 1 r3,8 = 0.0006 0.972 0.1

k4,8 = 0.2 P4,8 = 1 r4,8 = 0.0017

k5,8 = 0.2 P5,8 = 2.9 r5,8 = 0.0173

k6,8 = 0.2 P6,8 = 18.2 r6,8 = 0.0174

k7,8 = 0.2 P7,8 = 3.7 r7,8 = 0.0346

Детрит x(2) x9 7139.05 - - - - -

Биогены x(2) x10 280.35 - - - - -

Расчеты проведены в ППП Mathematica 5,0 [20], Задачи Коши решались методом Рунге-Кутта 4-го порядка. Для этого в среде Mathematica вызывается программный модуль

ClassicalRungeKutta /:

NDSolve'InitializeMethod[ClassicalRungeKutta,_]:=

ClassicalRungeKutta[];

А сам метод Рунге-Кутта реализован в виде процедуры:

ClassicalRungeKutta[_] ["Step" [f_, t_, h_, y_, yp_]]: = Block[{deltay, kl, k2, k3, k4>, kl = yp;

k2 = f[t + 1/2 h, у + 1/2 h kl] ; k3 = f[t + 1/2 h, у + 1/2 h k2]; k4 = f [t + h, у + h k3] ;

Рис, 2, Значения переменных и их фазовые портреты в первой модели (без учета промысла)

с1еКау = Ь (1/6 к1 + 1/3 к2 + 1/3 кЗ + 1/6 к4) ; {11, с1еКау}

];

где ур= Т (ж(га)), с1е^ау= ж(га+1) — ж(га), Ь - шаг. Для упрощения метода в процедуре используется фиксированный размер шага. Величина шага задается с помощью процедуры StartingStepSize, При численном моделировании в данной работе шаг 0.01

Рис, 3, Значения переменных и их фазовые портреты во второй модели (с учетом промысла)

Заключение

Описаны две модели потоков вещества в экосистеме Охотского моря. Системы дифференциальных уравнений построены па основе данных "неформальной"модели из |17| и уточнены по материалам |3|, На основе этой информации проведена идентификация параметров дня различных ситуаций, имитирующих отсутствие промысла и изъятие промысловых объектов. Показано, что точки равновесия систем уравнений пе являются устойчивыми. Проведено численное моделирование динамических процессов, возникающих в системе при отклонении от равновесного решения. Полученные результаты являются косвенным подтверждением того, что представления об экосистеме Охотского моря из |17| требуют уточнения. Были предприняты попытки "доработки", которые включали в себя введение дополнительного трофического блока биогенов, подбор констант ро и до, корректировку биомассы детрита (т^) и биогенов (т^), к = 1, 2, Однако это не привело к получению состояния равновесия систем уравнений пи дня одной из моделей.

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

Список литературы

[1] Горстко А.Б. Имитационная система "Азовское море". - Вопросы математического исследования и моделирования экосистемы Азовского моря, - Труды ВНИРО, 1976, т. CXVIII, с. 48-56.

[2] Домбровский Ю.А., Ильичев В.Г., Селютин В.В., Сурков Ф.А. Теоретические и прикладные аспекты моделирования первичной продуктивности водоемов. - Ростов-на-Дону: Изд-во РГУ, 1990,

[3] Дулепова Е.П. Сравнительная биопродуктивность макроэкосистем Дальневосточных морей. - Владивосток: Изд-во ТИНРО, 2001,

[4] Ивлев B.C. Экспериментальная, экология, питания рыб. - Киев: Наукова думка, 1977.

[5] Кучер А.И., Абакумов А.И. Рыбопродуктивность и динамика биомассы ихтиоценоза оз. Ханка. - Вопросы ихтиологии, 1997, т. 37, © 5, с. 619-624.

[6] Ладожское озеро - критерии состояния экосистем,ы. - СПб.: Наука, 1992.

[7] Леонов А.И., Осташенко М.М., Лаптева Е.И. Математическое моделирование трансформации органического вещества и соединений фосфора, в водной

среде: предварительный анализ условий функционирования экосистем,ы, Ладож-

©

[8] Ляпунов A.A. Об изучении балансовых соотношений в биогеоценозе (попытка м,а,тематического анализа). - Журнал общей биологии, 1968, т. 296, вып. 6, с. 25-32.

[9] Меншуткин В.В. Математическое моделирование популяций и сообществ водных животных. - Л.: Наука, 1971.

[10] Меншуткин В.В., Воробьев О.Н. Модель экологической системы Ладожского озера. - Современное состояние экосистемы Ладожского озера. Л.: Наука, 1987, с. 187-204.

[11] Понтрягин A.C. Обыкновенные дифференциальные уравнения. - М,: Наука, 1978.

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

[13] Сорокин Ю.И. Продукция миклофлоры. - Биологическая продуктивность океана. - М.: Наука, 1977, т.2, с. 209-233.

[14] Сорокин Ю.И. Черное море. Природа и ресурсы. - М,: Наука, 1982, 216 с.

[15] Сорокин Ю.И., Сорокин П.Ю., Сорокина О.В., Мамаева Т.И. Первичная продукция и гетеротрофный микропланктон в Охотском море. - Журн, общ, биологии, 1995, т. 56, Ж0- 5, с, 603-628,

[16] Сурков Ф.А. О динамике водных масс Азовского моря. - Вопросы математического исследования и моделирования экосистемы Азовского моря, - Труды ВНИРО, 1976, т. CXVIII, с. 56-62.

[17] Шунтов В.П., Дулепова Е.П. Современный статус, био- и рыбопродуктивность экосистемы, Охотского моря. - М,: Изд-во ВНИРО, 1997, с. 248-261.

[18] Эльсгольц Л.Э. Дифференциальные уравнения и вариационное исчисление. -М.: Наука. 1969.

[19] Schnute J.Т., Richards L.J. Use and abuse of fishery models. - Canadian Journal of Fisheries and Aquatic Sciences. Vol. 58, Ж0- 1, January 2001, pp. 10-17.

[20] Wolfram S. The Mathematica Book, J^th ed. - Wolfram Media, Cambridge University Press, 1999.