Алексеев В.В. Физическое и математическое моделирование экосистем

Подождите немного. Документ загружается.

Идентификация

экологических моделей имеет ряд существен-

ных особенностей. С математической точки зрения модели эколо-

гических процессов часто носят неустойчивый характер и с тру-

дом поддаются идентификации. Экспериментальная биологическая

информация

имеет значительные погрешности, намного превосхо-

дящие уровень шумов в физических системах, что также

услож-

няет

подбор коэффициентов модели. Неустойчивость системы

уравнений и большие ошибки экспериментальной информации за-

трудняют

декомпозицию модели, а если она

даже

произведена

каким-то

образом, то

результаты

данного этапа идентификации

часто не

подходят

для идентификации экосистемы в целом, т. е.

необходимо, как правило, настраивать параметры для всей мо-

дели в целом. При моделировании пространственных экосистем

заметные ошибки при определении коэффициентов

могут

вно-

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

Определение экологических параметров в значительной мере

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

уравнений, характера и объема экспериментальных данных. При

идентификации

имитационной модели необходимо в возможно

большей степени уменьшить неопределенность связей

между

па-

раметрами, исходя из физических и биологических соображений,

только после этого оставшиеся свободными параметры иденти-

фицировать

по данным наблюдений [50]. Важное значение для

уменьшения числа степеней свободы при выборе параметров

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

устойчивости особых точек и времени наступления пиков числен-

ности

в сезонной динамике компонент. Целесообразно начинать

процедуру

построения модели с качественной имитационной мо-

дели минимальной сложности, поведение которой на качествен-

ном

уровне описывает исследуемый эффект. После этого можно

производить дальнейшую детализацию модели, для достижения

необходимой точности количественного описания эффекта. Ис-

пользование предлагаемых процедур позволяет во многих слу-

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

Методы

идентификации,

использующие

решения

«прямых»

за-

дач. Данные методы идентификации основаны на определении

численных значений параметров модели, минимизирующих (на-

пример,

в среднеквадратическом смысле), отклонение

между

на-

блюдаемыми данными и расчетными значениями динамики си-

стемы.

Рассмотрим сначала простейший

случай

определения

коэффи-

циентов

логистического уравнения, описывающего динамику рас-

тущей

популяции [50, 103]:

-g- =

eX-

(10.26)

где X — плотность популяции; е, у — постоянные коэффициенты.

Решение

(10.26)

имеет вид

)e-

'-'

(10.27)

34!

где Х = е/у — параметр

«емкости»

среды; X(to)—значение чис-

ленности

популяции в начальный момент времени to. При малых

значениях X (на начальной стадии) рост численности микроорга-

низмов

носит экспоненциальный характер. Однако в дальнейшем

происходит постепенное замедление скорости роста и численность

организмов достигает предельного значения Х = г/у. Решая

(10.27)

относительно ехр [s(t— to)], получаем

(t-«.

)s=

X/X

(/„)-!

Х/Х(/)1

Логарифмируя последнее выражение, находим е = 1пФ(Х)/(£— to).

Используя метод наименьших квадратов для оценки

коэффици-

ента роста, окончательно получаем выражение

/(*

-*о).

(Ю.29)

где Xft — значение численности Х(^) в момент времени t — tu,.

п — число наблюдений.

Рассмотрим вопрос об определении коэффициентов классиче-

ской

системы «хищник-жертва»

- = е,Х, —

(10.30):

£X YXX

где Xi — плотность популяции жертвы; Хг — плотность популяции

хищника,

ei; е

; у и; Y21 — коэффициенты системы. Как показано

В. Вольтерра [45, 323], численности хищников и жертв в такой

системе совершают незатухающие колебания, причем изменения

численности хищника отстают по фазе от колебаний численности

жертв.

Пусть из данных наблюдений известны следующие величины

период колебаний численности компонент системы Т, отношение

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

ленности

жертв и хищников равны Х

и Хг. Тогда для опреде-

ления

четырех коэффициентов (ei, ег, Y12. Y21) имеем систему четы-

рех алгебраических уравнений:

Ц\

V ei

2 •

V12

Y21

_ (10.31):

Xi,

Х2.

342

Описанный

способ определения коэффициентов пригоден для си-

стемы хищник—жертва, испытывающей колебания вблизи поло-

жения равновесия.

В общем

случае

для идентификации коэффициентов может

быть использован метод, изложенный в

работах

[19, 53]. Проил-

люстрируем его на примере модели динамики многовидовой эко-

системы:

™

i/XiX,,

(10.32)

где X,- — биомасса /-го вида; е, и уц— коэффициенты модели. За-

пишем уравнение

(10.32)

в интегральной форме:

'* 'ft я 'ft

=e/

X,dt+

i/ i

X^X/d/, (10.33)

где ^o—начальное время наблюдений; th— время взятия пробы.

Поскольку

Xi(to)

и Xi(tk) известны из наблюдений, интегралы

(10.33)

могут

быть вычислены при помощи методов численного

интегрирования. Записывая уравнение

(10.33)

для различных ин-

тервалов времени (количество которых равно числу неизвестных

коэффициентов

уравнения), получаем систему линейных алгеб-

раических уравнений, решая которую можно найти е* и Y;J-

Вследствие

большого разброса экспериментальных значений, ко-

эффициенты

модели

будут

изменяться в зависимости от того, ка-

кие

интервалы (to, th) взяты для их вычисления. Чтобы сгладить

эти

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

коэф-

фициентов,

составим квадратичную форму:

Г h *к п

k "I

Ф(е„

уц) = £ Р* X, - г, \ X, di + £

</ \ X,X/ dt . (10.34)

Числа рй выбираются в зависимости от дисперсии ошибок изме-

ряемых величин, значениям с большей дисперсией ошибки ста-

вятся в соответствие меньшие значения р&. Поскольку при изу-

чении динамики системы наиболее важно определить положение

максимальных величин компонент, соответствующие значения р^

следует

брать наибольшими. Для натурных значений (данных),

полученных через одинаковые промежутки времени с примерно

одинаковой дисперсией ошибки, в качестве коэффициентов р& мо-

жно взять отношения измеряемых в момент th значений числен-

ности к их среднему значению. Приравнивая к

нулю

первые про-

изводные положительно определенной квадратичной формы

получаем систему линейных алгебраических уравнений, решая ко-

торую,

находим значения zi и y

, минимизирующие Ф.

343

При

идентификации сложных имитационных моделей опреде-

ление параметров обычно производится

в ряд

этапов.

На

первом

этапе

по

возможности, значения эколого-физиологических

коэф-

фициентов определяются экспериментально. Однако необходимо

учитывать,

что

коэффициенты, определенные

условиях лабора-

торных

(и

даже

натурных) экспериментов

отдельными компо-

нентами экосистемы

могут

отличаться

от

значений параметров

в природных условиях

для

экосистемы

целом.

дальнейшем,

полученные экспериментально значения коэффициентов подстав-

ляются

уравнения модели

качестве начального приближения

при

многократных просчетах

на ЭВМ

(имитациях) уточняются

до

тех пор,

пока

не

будет

достигнуто приемлемое

(в

пределах

точ-

ности определения натурных данных) соответствие

между

экспе-

риментальными

расчетными значениями. Например,

работе

[244]

по

описанной

схеме

были определены параметры модели

эвтрофикации

озера Онтарио, идентифицированной

по

данным

наблюдений

1967—1970

гг. Для

достоверной идентификации

мо-

дели было сделано

вариантов расчета. Показано,

что

модель

хорошо воспроизводит весенний

пик

численности фитопланктона,

последующее

его

отмирание

середине лета

осеннее цветение.

Метод

интегральных преобразований. В. Б.

Георгиевским

[50]

был

предложен способ идентификации экологических моде-

лей, основанный

на

использовании метода интегральных преобра-

зований.

При

использовании этого метода производится исключе-

ние

производных

результате

замены дифференциального урав-

нения

некоторым интегральным алгоритмом, приводящим

к си-

стеме алгебраических уравнений относительно идентифицируемых

параметров.

частности, одно

из

преобразований, позволяющее

исключить

из

модели

L = 0

производные первого порядка, имеет

вид

\ \Ldt = 0.

(10.36)

Для определения параметров

ei, ег, Y12, Y21

системы хищник—

жертва

помощью

(10.36)

можно составить систему алгебраиче-

ских уравнений:

[££] е,С/

[X,]

1Я

[Х,Х

[-££-]

= -г

[X,]

Yll

£/ [Х.Х.],

где интегральное преобразование осуществляется

по

четырем

об-

ластям.

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

ляет заменить некорректную операцию дифференцирования экспе-

риментальных функций операцией интегрирования,

при

которой

случайные экспериментальные ошибки сглаживаются.

Эта

мето-

дика идентификации применима

как для

обыкновенных дифферен-

344

циальных уравнений,

так и для

уравнений

частных производ-

ных. Весьма существенно,

что

идентификация параметров

не тре-

бует

предварительного решения краевых задач.

Анализ

чувствительности.

Анализ экспериментальных данных

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

могут

производиться

при

помощи мето-

дов теории чувствительности

[103, 318],

которая позволяет усовер-

шенствовать процедуру перебора параметров

ввести

рассмот-

рение

некоторые качественные характеристики экосистем, мало

зависящие

от

численных значений параметров. Теория чувстви-

тельности позволяет интерпретировать экологические процессы

с точки зрения выявления тенденций

в их

поведении.

Пусть имеется система уравнений

(10.37)

X(t

)=a,

где

X —

вектор состояния экосистемы;

р —

вектор параметров;

X(to)—вектор начальных состояний.

Под

функциями чувстви-

тельности понимаются производные

от

компонент вектора состоя-

ния

по

параметрам. Дифференцируя

(10.37)

по

параметрам

начальным условиям, получим

(дх,\_

df,

dXk

dt \ dpi )~ dX

dpi "т"

(10.38)

dpi

•

Аналогично

могут

быть получены выражения

для

функций

чувст-

вительности более высокого порядка. Например,

для

функций

чув-

ствительности второго порядка

по

параметрам

имеем

_ dfi

. d'fi dX

d*fi

dp,-

dpi dXk dpi dpi '

dXk

dpi dpj "• dpj dp

(10.39)

dpj

dpi

Методы теории чувствительности

могут

быть использованы

для распределенных

пространстве моделей экосистем.

позиций теории чувствительности экологической моделью

является система уравнений

(10.37)

вместе

уравнениями

для

•функций

чувствительности (10.38), (10.39). Поведение функций

345

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

«узкие»

звенья системы,,

определить пространственно-временные интервалы наибольшей

уязвимости (по экстремуму функции), а также наличие или от-

сутствие взаимодействия (синергизма) факторов, управляющих

поведением системы. Во многих случаях функции чувствительно-

сти

могут

использоваться в качестве наглядных интегральных ха-

рактеристик возможных тенденций в динамике экосистемы. Функ-

ции

чувствительности первого порядка показывают тенденциг

росту или уменьшению экологических переменных (численности

или

биомассы) при изменении значения исследуемого параметра..

Влияние параметра является максимальным в области экстремума

Рис.

10.14.

Функции чувствительности первого (а) и второго (б) поряд-

ков

для логистической модели роста (10.26).

функции

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

систему

отсутствует

в области, где функция чувствительности

равна нулю. Функции чувствительности по начальным состояниям

позволяют выявить эффекты последействия, а также моменты наи-

большего влияния изменений начального состояния на динамику

системы. Смешанные производные — функции чувствительности

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

параметров (в том числе и антропогенных факторов). Заметим,

что анализ чувствительности модели конструктивен при сравни-

тельно небольшом числе экологических параметров. В качестве

примера на рис. 10.14 представлены функции чувствительности

для логистической модели (10.26). Как видно из рисунка, макси-

мальная чувствительность при воздействии на параметр роста е

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

ный

момент времени. При исследовании изменений параметра

•у достаточно исследовать чувствительность в стационарной точке

346

Х = е/у. Анализ функции чувствительности второго порядка пока-

зывает, что предыстория системы (чувствительность по началь-

ному состоянию а) существенна при воздействии на оба пара-

метра логистической модели. В стационарной области зависимость

от начальных условий исчезает. Имеет место синергизм воздейст-

вий

на параметры логистической модели е и Y-

Следует

заметить,

что при нормировании антропогенного воздействия на скорость

роста в с

учетом

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

метра Де заметно меньше изменений этого параметра, оценивае-

мых в стационарной области.

Идентификация

параметров

модели

экосистемы

водоема-

охла-

дителя.

Идентификация параметров модели экосистемы водоема-

охладителя производилась на ЭВМ на основе метода чувствитель-

ности,

при помощи разработанной В. Н. Михайловым универсаль-

ной

диалоговой системы «Идентификатор» [83]. Пусть динамика

системы описывается в общем виде векторным дифференциальным

уравнением

JX-

F{t, X, W, p), X(t

) = X

, (10.40)

-+•

-*•

где X — n-мерный вектор состояния системы; W — ^-мерный век-

тор входных возмущений; р — r-мерный вектор параметров;

F(t, X, W, р) —нелинейная вектор-функция. Пусть имеется также

вектор измерений У; (размерности р^я). Требуется по измере-

ниям

Yi(t) при 0 =SS t ^ T найти вектор определяемых парамет-

ров р. Алгоритм идентификации состоит в следующем. Вектор

переменных состояния системы представляется в виде

X(t, p) = X(t, P*) + v(t,

p*)bp,

(10.41)

где р* — вектор начальных значений параметров; v{t, р*)—мат-

рица

размерности (пХг), элементами которой являются функции

чувствительности вектора X по параметрам р^

(н-=1>

О- Тогда,

считая матрицу v(t, p*) известной, вектор р находим итерацион-

ной

процедурой по формуле

р = р* + Х/±р,

(10.42)

где Ар в каждой итерации определяется в

результате

наилучшего

приближения

к измеряемому вектору У; по

методу

наименьших

квадратов. Процедура

(10.42)

продолжается до получения резуль-

тата с требуемой точностью. Оценку погрешности можно произ-

водить, например, по разности

между

векторами настраиваемых

параметров или значениями функционала невязки выходных экс-

периментальных и моделируемых переменных, полученных в

двух

следующих

друг

за

другом

итерациях.

Идентификация

параметров экосистемы водоема-охладителя

проводилась для упрощенной модели, включающей две экологи-

ческие группы микроводорослей (ASP и ASZ), зоопланктон (М4),

-концентрации

в воде лимитирующего биогенного элемента (фос-

347

фат

Sp) и

детрита

(M<j) [см.

уравнения (10.15), (10.16), (10.18),

(10.19), (10.20)].

Для

идентификации модели использовалась

сле-

дующая экспериментальная информация

—

данные наблюдений

за численностью отдельных видов фитопланктона, зоопланктона,

концентрацией фосфатов

воде,

также

за

входными парамет-

рами системы (температурой

световой экспозицией). Уравнения

модели интегрировались

по

времени методом Рунге—Кутта

с по-

стоянным шагом, равным

0,2 сут.

Начальные значения идентифи-

цируемых параметров приведены

табл.

10.4. Там же

даны инте-

Таблица

10.4

Начальные

оценки параметров

интегральные

значения

функций

чувствительности модели экосистемы водоема-охладителя

Параметр

eoi

Ppi

= Yi

Pel

Yd4

Уй

602

Y24

= Y2

Рс2

Р.

Р>

V-o

£04

6з

Начальная

оценка

0,26

0,00273

0,01

0,0452

0,01

0,352

0,008

0,055

0,6

0,00008

0,2

0,15

0,08

0,045

"ASP

5,2-Ю

7,6-10

8,0-Ю

3,7-10»

0,066

3,8-Ю

1,4.10

1,4-Ю

4,9-10

8,0-!0

1,4.10

0,024

0,3

1,9-10

Функции

чувствительности

"ASZ

5,4-Ю

1,4-10»

2,9-Ю

3,6-10»

6,7-10

4,7-Ю

5,3-Ю

2,7-Ю

3,8-10

7,5-10

8,0-10"

7,9-Ю

1,9-Ю

2,9-Ю

6,4-10'

6,6-10

0,5

2,9-Ю

0,013

1,3-Ю

5,0

8,5-10

0,18

3,9

2,2-10»

2,4-Ю

9,6-Ю"

0,022

0,24

"зоо

1,6-Ю

1,2-10

4,0-Ю

9,5-10

2,0-Ю

3,4-Ю

9,9-10

2,5-Ю

0,92

1,3-10

9,7-10

3,1-10

6,2-10

тральные значения функций чувствительности

по

каждому изме-

ряемому

выходу,

которые определялись следующим образом:

dt,

(10.43)

(»'=

1, n; 1=1, r),

где y'i =yi(pi, Pi, .... pi-i, pi + bpi, pi+i, ..., pr); 6pi=0,0\ p

;

—

порядковый номер измеряемой переменной (выхода)

си-

стемы; /

—

номер определяемого параметра;

п —

число перемен-

ных;

г —

общее число параметров; 0,-

—

весовые коэффициенты,

значения

которых обратно пропорциональны среднеквадратиче-

ским

ошибкам измеряемых переменных. Значения некоторых

па-

раметров модели фиксировались

до

идентификации, исходя

из фи~

зических соображений.

частности, полагались фиксированными-

значения

показателей экспонент процессов биосинтеза

распада

348

в температурных зависимостях скоростей роста организмов (а

0,065

(С

«р = 0,2

(С

см. п. 2.7), что значительно упро-

щало

процедуру

идентификации, которая становилась линейной

задачей по отношению к идентифицируемым параметрам.

Процедура идентификации модели была разбита на два этапа.

Сначала идентифицировались параметры, входящие в уравнения

динамики численности фитопланктона, концентрации минераль-

ного фосфора и

детрита.

При этом в качестве

входных

перемен-

ных задавались биомасса зоопланктона, температура воды и све-

товая экспозиция, а выходными измеряемыми переменными яв-

лялись численности микроводорослей фитопланктона и концентра-

ция

фосфора. На втором этапе производилась идентификация

параметров, входящих в уравнение динамики зоопланктона. В ка-

честве

входов

использовались значения температуры и световой

экспозиции;

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

биомасса зоопланктона; значения параметров, определенных на

первом этапе, задавались фиксированными. Окончательная на-

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

одновременно с использованием

всех

выходных

измеряемых пере-

менных. При этом малоинформативные (низкочувствительные) па-

раметры

могут

быть зафиксированы. Анализ чувствительности

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

являются параметры температурных зависимостей скоростей ро-

ста фитопланктона обеих групп микроводорослей (p

i, Pc2, {5

8oi, еог), а также коэффициент высвобождения фосфора при раз-

ложении

детрита.

Наиболее чувствительными переменными модели

являются численности водорослей планктона, причем численность

микроводорослей летне-осеннего комплекса наиболее сильно зави-

сит от интенсивности минерализации

детрита.

Фитопланктон яв-

ляется также более чувствительным к изменению коэффициентов

фотосинтеза (p

i, Рсг), чем смертности (e

i и еог). Изменения био-

массы зоопланктона наиболее чувствительны к параметрам роста

фитопланктона, коэффициентам потребления микроводорослей и

смертности зоопланктона. Содержание в

воде

лимитирующего био-

генного элемента (фосфора) наиболее чувствительно к изменениям

параметров фотосинтеза фитопланктона

минерализации дет-

рита. Параметры потребления фитопланктона и

детрита

зоопланк-

тоном (уь Y2. yd), содержания фосфора в химическом составе мик-

роводорослей (р

), осаждения

детрита

), усвоения зоопланк-

тоном микроводорослей (U

) в основном являются

малочувст-

вительными. Выполненный аналз позволяет считать, что пара-

метры роста фитопланктона, по-видимому,

могут

быть выбраны

в качестве наиболее чувствительного экологического индикатора

по

отношению к воздействию антропогенных факторов.

СПИСОК

ЛИТЕРАТУРЫ

1. Абросов Н. С, Боголюбов А. Г. Экологические и генетические за-

кономерности

сосуществования и коэволюции видов.— Новосибирск: Наука.

1988.—273

с.

2. А б р о с о в Н. С, Ковров Б. Г. Анализ видовой структуры трофиче-

ского уровня одноклеточных.— Новосибирск: Наука.

1977.—

190 с.

3. А б р о с о в Н. С, Ковров Б. Г., Черепанов О. А. Экологические

механизмы сосуществования и видовой регуляции.— Новосибирск: Наука.

1982.—301

с.

4. Абросов Н. С. Замкнутая экосистема и

хемостат

— аналогия и примеры

анализа//Журн. общ. биол.—

1979.—

Т. 40, вып. 2.— С.

252—258.

5. Айзатуллин Т. А., Шамардина И. П. Математическое моделиро-

вание экосистем континентальных водотоков и водоемов//Итоги науки и

техники. Общая экология. Биоценология. Гидробиология.— М.: изд.

ВИНИТИ,

1980.—С.

154—228.

6. Алек и н О. А., Л я хин Ю. И. Химия океана.— Л.: Гидрометеоиздат,

1984.—343

с.

7. Алексеев В. В. О роли консументов в формировании разнообразия

видов в биогеоценозах//Журн. общ. биол.— 1973. Т. 34, вып. 3.— С. 451—

458.

8. А л е к с е е в В. В. Влияние фактора насыщения на динамику системы

«хищник—жертва»//Биофизика,—1973. Т. 18. № 5,—С.

922—926.

9. Алексеев В. В. Исследования модели двухкомпонентных фитоценозов//

Журн. общ. биол.— 1975. Т. 36, вып. 6.—С.

914—919.

10. Алексеев В. В. Биофизика сообществ живых организмов//Успехи физи-

ческих наук—1976. Т. 120. Вып. 4.—С.

647—676.

11.

Алексеев В. В. Динамические модели водных биогеоценозов//Человек

биосфера. Вып. 1—М.: МГУ, 1976,—С.

3—137.

12. Алексеев В. В. Биогеоценозы-автогенераторы и триггеры//Журн. общ.

биол.— 1977. Т. 37, Вып. 6.—С.

738—744.

13.

Алексеев В. В. Кинетика неравновесных процессов в экологии//Термо-

динамика

и кинетика биологических процессов.— М.: Наука,

1980.—

С. 235—

240.

14. Алексеев В. В. Теоретические и экспериментальные исследования эколо-

гических систем//Итоги науки и техники. Общая экология. Биоценология.

Гидробиология. Т. 5. Моделирование водных экосистем.— М.: изд.

ВИНИТИ,

1980.—С.

58—104.

15. Алексеев В. В. Изменение числа видов в экосистеме в зависимости

от температуры и возможные последствия теплового загрязнения в связи

с развитием энергетики//Экологические аспекты исследований водоемов-

охладителей АЭС. М.: изд. ИЭМЭЖ АН

СССР,

1983.—С. 170—181.

16. Алексеев В. В., Кокорин А. О. Математическая модель влияния

пестицида на биоценоз//Тр. ИЭМ. Вып. 12. (98).—1982.—С.

111—116.

17. Алексеев В. В., Корниловский А. Н. Автостохастические процессы

биофизических системах//Биофизика.—1982. Т. 27. Вып. 5.—С. 890.

18. Алексеев В. В., Крыш ев И. И. Кинетические уравнения для описа-

ния

динамики биоценозов//Биофизика.— 1974. Т. 19, вып. 4.— С.

754—759.

350