Discrete-continuous model for the problem of analysis critical level of exploitation of bioresources

Full Text

Введение. Традиционно математические методы в промысловой ихтиологии были связаны со статистическим анализом. Проводился поиск корреляционных взаимосвязей по данным об изменениях факторов среды и ряда величин, характеризующих популяционные процессы. Составлялись уравнения регрессии, прогностическая значимость которых нивелировалась тем, что значения популяционных характеристик редко удаётся определить с достаточной точностью или они могут выражаться только в относительных показателях. Работа Роберта Мэя [1] представила новое направление для подобных исследований и привлекла внимание биологов к возможностям применения нелинейных динамических моделей, исследуемых с применением компьютерных сред. Практические задачи моделирования в данной области имеют такую особенность, что при обосновании адекватности моделей необходимо дать непротиворечивую интерпретацию изменений в их качественном поведении. Разнообразие нелинейных эффектов предполагает выявление возможных результатов, логически отбрасываемых при переходе обратно к биологической терминологии. Для каждой моделируемой ситуации необходимо определить набор выводов, которые допустимо сформировать в итоге вычислительного исследования. ISSN: 2310-7081 (online), 1991-8615 (print); doi: http://dx.doi.org/10.14498/vsgtu1116 © 2014 Самарский государственный технический университет. Образец цитирования: А. Ю. П е р е в а р ю х а, “Дискретно-непрерывная модель в задаче анализа критического уровня эксплуатации биоресурсов” // Вестн. Сам. гос. техн. унта. Сер. Физ.-мат. науки, 2014. № 1 (34). С. 145-155. 145 А. Ю. П е р е в а р ю х а В настоящей статье рассматриваются динамические системы, моделирующие важные изменения популяционных процессов. Обсуждается вопрос о том, насколько свойства нелинейных функциональных итераций соответствуют задачам и представлениям популяционной проблематики. Решается обусловленная анализом ретроспективных данных о величине уловов осетровых задача математической формализации эмпирически выявленной закономерности в виде функциональной зависимости численности нерестового запаса и эффективности воспроизводства. Обосновывается сценарный подход к моделированию последствий эксплуатации биоресурсов. 1. Дискретные модели теории воспроизводства популяций. Для моделирования в ихтиологических исследованиях фундаментом стала теория формирования пополнения популяций, предложенная У. Рикером в [2] в виде опирающейся на реальные данные систематической концепции. Анализ результатов предпринятого Норвежским институтом морских исследований комплексного моделирования экосистемы Баренцева моря показывает важную роль в большой модели блока воспроизводства. Проблема математической формализации в рамках теории о регулируемой зависимости между нерестовым запасом и формирующимся пополнением актуальна для исследования популяционных изменений рыб, мигрирующих из морей в реки для размножения. Тогда наблюдатели могут точнее подсчитать численность родительского запаса. Можно определить площадь постоянно пригодных нерестилищ при известном уровне реки (для зарегулированной Волги сток и температура в половодье измеряются точно) и при систематических мониторинговых наблюдениях за скатом молоди в створе реки оценить численность нового поколения. У. Рикер [2] предложил часто применяемую функцию R = f (S) связи «запас→пополнение», отражающую негативное действие на выживаемость поколения повышенной плотности запаса. Будем рассматривать эту модель как дискретную динамическую систему в виде итераций {ψ (j) }j 0 , где R0 , R1 , R2 , . . . - последовательность точек, описывающих эволюцию системы и определённых условием Rj+1 = aRj exp(-bRj ). Её свойства как функции и недостатки как популяционной модели обсуждались неоднократно, например в [3]. Биологический смысл модели становится яснее, если перейти к дифференциальному уравнению убыли текущей численности N (t): dN = -(αN (0) + β)N (t), dt t ∈ [0, ϑ]. (1) Константы заданного на промежутке времени обыкновенного дифференциального уравнения (ОДУ) соотносятся с константами формулы Рикера: a = λ exp(-βϑ), b = αϑ, где λ - средняя плодовитость популяции. Согласно теории, повышенная плотность на нерестилищах становится негативным фактором, увеличивающим смертность. Предполагается, что на смертность влияет именно начальная 146 Дискретно-непрерывная модель в задаче анализа критического уровня эксплуатации биоресурсов численность поколения. Большое значение имеет длина интервала ϑ - периода уязвимости поколения по отношению к основным факторам смертности. Мы предлагаем на практике для вариативности в оценках параметров использовать модель именно в виде численного решения ОДУ (1) с начальным условием N (0) = λS. Динамическая система ∀ j > 0 обладает аттрактором A таким, что ∀R0 lim {ψ (j) } = A. j→∞ Для простого случая аттрактор глобальный и является устойчивым состоянием равновесия: A ≡ R∗ = ψ(R∗ ). Первый метаморфоз A определён условием нарушения критерия устойчивости неподвижной точки |ψ (R∗ )| < 1, следующего из теоремы, независимо доказанной Д. Гробманoм и Ф. Хартманoм в [4, 5]. При значении производной ψ (R∗ ) = -1, когда a = e2 , происходит бифуркация удвоения периода. У функциональной итерации возникают две новые циклические точки ψ n (R∗ ) = ψ n+2 (R∗ ), являющиеся неподвижными точками второй итерации ψ 2 (R). При дальнейшем увеличении параметра две неподвижные точки ψ 2 (R) аналогично потеряют устойчивость и появятся n = 4 точек ψ 4 (R). Любая ранее устойчивая R∗ становится репеллером и делит область притяжения образовавшегося цикла на субинтервалы Υ1 , Υ2 . Для R0 ≈ R0 , R0 ∈ Υ1 , R0 ∈ Υ2 , отличается порядок обхода цикла ∗ lim ψ(R0 )n = R1 при n = 2i n→∞ и, соответственно, ∗ lim ψ(R0 )n = R1 при n = 2i . n→∞ При изменении управляющего параметра в диапазоне значений e2 < a < a ˆ реализуется каскад бифуркаций удвоения периода, происходящих при ψ 2n (R∗ ) = -1. Когда период цикла становится бесконечным, образуется странный аттрактор, гомеоморфный канторовскому множеству. Скорость увеличения периода цикла описывается универсальными константами М. Фeйгeнбаума [6], в честь которого называют данный сценарий хаотизации. Из-за бесконечного количества малых субинтервалов наблюдается чувствительная зависимость от начальных условий. Практические цели предметной области диктуют необходимость биологической интерпретации сценария хаотизации. Бифуркационный параметр a в модели Рикера трактуется как мера репродуктивного потенциала популяции, что приведёт исследователя к очевидным заключениям. Мы отметим ряд проблем для поиска практических решений при моделировании. Во-первых, сценарий хаотизации связан с возникновением 147 А. Ю. П е р е в а р ю х а ряда других нелинейных эффектов. Метаморфозы происходят как с аттрактором, так и с границами его области притяжения, и подобным явлениям трудно дать популяционную интерпретацию. Во-вторых, существуют альтернативные модели вида «запас - пополнение», в которых хаотизация происходит при увеличении параметров, трактуемых иначе. В-третьих, наблюдается огромная, маловероятная для реальных популяционных процессов, амплитуда хаотических колебаний. Сценарий хаотизации является свойством определенного класса дискретных динамических систем SU -отображений, значит, необходимо вести речь о биологической интерпретации свойств отображений, удовлетворяющих критериям теоремы Д. Сингера [7]. Задачи промысловой ихтиологии пересекаются с фундаментальными и ещё не полностью разрешёнными проблемами нелинейных эффектов, т. к. в реальных моделях может присутствовать несколько бифуркационных параметров. 2. Модель регулируемого плотностью воспроизводства. Характер зависимости решения (1) между запасом и пополнением с резко выраженным единственным экстремумом далеко не всегда согласуется с данными наблюдений ихтиологов и с точки зрения математики означает возможность появления циклов всевозможных периодов в соответствии с порядком теоремы А. Н. Шарковского [8, 9]. Опубликованные результаты наблюдений [10] свидетельствуют о возможности появления двух максимумов Sm1 , Sm2 на графике зависимости нерестового запаса и пополнения. Очевидный и простой способ описания сложной формы кривой - воспользоваться второй итерацией формулы Рикера f (f (S)) ≡ f 2 (S). Подобные попытки моделирования известны. Однако трудно представить, что выполнится f (Sm1 ) = f (Sm2 ), более того, данный способ некорректен из-за свойств функции R = f (S): df (S) = ae-bS (1 - bS), dS df (S) = a; S→0 dS lim df 2 (S) = a2 exp (-bS - abSe-bS )(1 - abSe-bS - bS + ab2 S 2 e-bS ), dS df 2 (S) lim = a2 , S→0 dS подразумевающих, что при деградации запаса скорость воспроизводства стремится к предельному значению a, но в f 2 возрастает уже до a2 . По данным о нересте осетровых рыб неоднократно отмечалось, что при существенном снижении плотности запаса эффективность воспроизводства резко уменьшается. Иначе происходят бифуркации Rn+1 = f (f (Rn )), f (R) = aRe-bR . После потери устойчивости возникает не цикл, а аттракторы A1 , A2 , с которыми происходят удвоения периода. Задачу описания разнообразных сложных форм зависимости мы предлагаем решать при помощи введения функционалов с подобранной ограниченной областью значений E(Θ) = [l, 1), l > 0 в правую часть (1). 148 Дискретно-непрерывная модель в задаче анализа критического уровня эксплуатации биоресурсов Анализ опубликованных данных о воспроизводстве находящейся на грани выживании популяции волжского осетра позволил сделать вывод о сильном действии эффекта Олли [11]. Это описанное в 30-годы XX века зоологами нелинейное явление выражается в том, что если плотность популяции в большом ареале обитания становится ниже оптимальной, то непропорционально резко падает эффективность воспроизводства. Оценка смещения максимума распределения веса молоди в прудах при увеличенной сверх норматива плотности позволил нам предложить модель в виде системы ОДУ с функционалом Θ(S), описывающую убыль поколения на интервале уязвимости t ∈ [0, T ]: dN = - αw(t)N (t) + Θ(S)β N (t), dt g dw 1 = k , , Θ(S) = dt 1 - exp(-cS) N +ζ (2) где S - величина нерестового запаса; w(t) - функция, отражающая условный уровень размерного развития поколения, влияющий на увеличение пищевых потребностей; g - параметр, учитывающий ограниченность количества доступных кормовых объектов; Θ(S) - некоторая убывающая функция (Θ(S) → 1 при S → ∞), не влияющая на вычисление N (T ); T = 45 сут., если численность запаса достаточно велика. Введение в систему ОДУ (2) быстро убывающей функции Θ(S) отражает снижение эффективности воспроизводства при деградации популяции, связанное с уменьшением вероятности встречи особей в местах размножения. Влияние данного эффекта можно оценить по имеющимся данным о состоянии осетровых [12]. Далее в (2): ζ - параметр, учитывающий ограничение темпов развития, не зависящее от численности популяции; c - параметр, характеризующий степень выраженности эффекта Олли; α - мгновенный коэффициент компенсационной смертности; β - мгновенный коэффициент декомпенсационной смертности. Начальные условия для системы (2) определяются следующим образом1 : √ (3) N (0) = λS, w(0) = w/ λ. ¯ Типичный график моделируемой зависимости запаса и пополнения R = = f (S) = N (T ) (кривая 1 на рис. 1), полученный при численном решении задачи Коши (2), (3) в инструментальной среде AnyLogic, является унимо∗ ∗ дальной кривой с двумя нетривиальными пересечения R1 , R2 с биссектрисой N (T ) = S (прямая 2 на рис. 1). Кривая 1 отличается пологой ниспадающей правой ветвью с ненулевой горизонтальной асимптотой и крутым наклоном восходящей ветви только после первого пересечения с биссектрисой коорди∗ натного угла R1 < R < Rmax . Система (2) форме дискретно-непрерывной системы может быть использована в вычислительных средах, которые поддерживают многошаговые численные методы и содержат средства имитационного моделирования. В таких средах мы можем получить последовательность непрерывных систем, связанных в момент переопределения начальных условий Ni+1 (0) = Ni (T )λ, w(0) = const, 1 Биологи заметили автору, что начальный размер икринки должен уменьшаться при б´льшей плодовитости. о 149 А. Ю. П е р е в а р ю х а Рис. 1. Типичный график R = f (S) (кривая 1), получаемый при численном решении задачи Коши (2), (3) для t = T , и биссектриса координатного угла (линия 2) [Figure 1. Typical graph of R = f (S) (curve number 1) obtained from the numerical solution of the Cauchy problem (2), (3) with t = T , and image of coordinate angle bisector (line number 2)] и задающих дискретную траекторию. При организации модельных экспериментов необходимо вводить соответствие условного модельного времени отрезков t ∈ {[0, T ]i } календарному. Использование событийно-управляемого подхода позволяет выделить некое особое, приводящие к переменам состояние, что в сочетании с непрерывным описанием популяционного процесса перспективно при рассмотрении последствий изменяющегося по некоторым условиям внешнего воздействия. 3. Сценарное моделирование управления промыслом. Особенности новой модели удобно использовать при оценках последствий различных антропогенных воздействий на благополучие популяций. Для получения сравнительных характеристик различных стратегий эксплуатации биоресурсов мы предлагаем исследовать в вычислительной среде парные наборы модельных сценариев вмешательства в саморегулируемые популяционные процессы. Качественный сценарий будет отражать алгоритм изменений управляющего воздействия, решение о котором может быть принято в реальности контролирующим уровень воздействия экспертом. Наиболее интересна последовательность решений и реакции на них, приводящая к деградации биоресурсов. Так, можно исследовать эволюцию состояния с момента, когда давление промысловой эксплуатации начинает превышать неистощительный уровень изъятия. Развитие моделируемой ситуации при изменении режима промысла зависит от свойств фазового портрета динамических систем с предложенным оператором эволюции. Поведение дискретно-непрерывной системы с оператором эволюции на основе численного решения (2) качественно отличается от итераций функций Рикера и Шепарда разделением фазового пространства на две области притяжения. Граница областей ∂Ω1 представляет простой случай для муль∗ тистабильных систем, являясь неустойчивой точкой равновесия ∂Ω1 ≡ R1 , или «репеллером». Анализ устойчивости неподвижных точек динамической системы, реализованной в инструментальной среде моделирования, можно проводить с использованием свойства второй итерации ψ 2 (x). Необходимым 150 Дискретно-непрерывная модель в задаче анализа критического уровня эксплуатации биоресурсов и достаточным условием устойчивости неподвижной точки x∗ одномерного отображения являются неравенства ψ 2 (x) > x при x < x∗ и ψ 2 (x) < x при x > x∗ . Устойчивым является тривиальное состояние равновесия, таким образом, модель подразумевает наличие критически низкой численности для существования биологического вида. Интересным представляется поведение тра∗ ектории в окрестности точки репеллера R1 ±ε, когда под внешним воздействием, характеризующим некоторой неопределённостью, траектория медленно покидает эту окрестность. Рассмотрим сценарий, когда после длительного периода получения стабильных уловов возникает предположение, что запас популяции Sst недоиспользуются промыслом. Пример подобных заключений экспертов о промысле осетровых в Волге - статья [13]. Отметим, что учёт запасов специалистами производится по количеству пропущенных на нерест особей, тогда как статистика промысла ведётся по массе уловов, что может вносить искажения в оценки имеющихся данных. Разработанная модель пополнения оперирует единицами численности. Пусть первоначально установленная доля изъятия ¯ H0 в момент t увеличивается: H1 = H0 + ∆H, ∆H ≈ 0.15. В дискретно-непрерывной системе начальные условия теперь сопрягаются так: Ni+1 t=0 = (1 - H1 )λNi t=T . i Несколько следующих сезонов уловы превышают среднемноголетние за предыдущий пятилетний период. Эксперты, обосновавшие возможность увеличения уровня эксплуатации, утверждаются в рациональности своего решения. В сцепленных сценариях будем моделировать популяции с различной введённой дополнительной промысловой нагрузкой и последующим её изменением в некоторый момент времени, получая на графике две временные диаграммы величины уловов. На диаграммы на ось t, где отложена шкала времени инструментальной среды, так же выводится пересчёт модельного времени в нумерованные календарные сезоны Yi (одному сезону соответствуют 8 тактов модельного времени AnyLogic 5). Для удобства наблюдения отсчёт ¯ запускается от момента t, когда был изменен режим промысла. Далее принятие решений и развитие процесса может происходить разными путями, отражёнными в предикативно заданных условиях двух модельных экспериментов с одинаковой исходной численностью запасов при H0 = 0.4, λ = 227000.2 После кратковременного увеличения уловы снижаются до прежних значений и продолжают медленно уменьшаться со скоростью, зависящей от ∆H. Эксперты предполагают стабилизацию популяции, не вносят изменения в режим промысла и далее происходит резкое падение уловов, после которого остановка лова не спасает от деградации (кривая 1 на рис. 2): Ri 2 Sst . Взята оценка среднемноголетней плодовитости каспийской севрюги. 151 А. Ю. П е р е в а р ю х а Рис. 2. Сценарии деградации при различном уровне перелова [Figure 2. Degradation scenarios at various level of overfishing] Длительный сценарий деградации растянулся на 33 модельных сезона. Предположим, для второй популяции H1 значительно превышает оптимальный уровень эксплуатации. Объёмы вылова падают быстрее, и эксперты, осмыслив негативные последствия, решают через 21 сезон вернуться к уровню H0 , ∗ когда состояние популяции находится в ε окрестности репеллера R1 (псевдостабильно). Некоторое время уловы постоянны, но незначительные изменения смертности неожиданно резко выводят популяцию из нестабильного равновесия, и происходит деградация (кривая 2 на рис. 2). Поведение популяции низкой численности перестаёт удовлетворительно описываться детерминированной динамикой средних величин, и мы получаем очень мало сведений о реальном состоянии потерявших промысловое значение биоресурсов. В случае, когда эксперты после получения данных о сокращении уловов за предыдущие 3 сезона своевременно переводят промысловый режим H2 = ∗ = H0 - ∆H при Rj > R1 , обе имитированные в сценариях популяции за 39 модельных сезонов восстанавливают свои промысловые запасы (рис. 3). Длительный сценарий деградации (подобный кривой 1 на рис. 2) реализо- Рис. 3. Сценарии восстановления после прекращения перелова [Figure 3. Restore scenarios because of a stop overfishing] 152 Дискретно-непрерывная модель в задаче анализа критического уровня эксплуатации биоресурсов вался для популяций осетровых Каспийского моря. В 1979 г. уловы достигли рекордных 27 тыс. тонн. В 1988-89 гг. произошло резкое падение уловов, промысел не был приостановлен вовремя, эксперты рассчитывали на успех технологии искусственного воспроизводства. Только в ноябре 2010 г. по инициативе Казахстана главы пяти государств Каспийского бассейна одобрили введение пятилетнего моратория на промышленный лов осетровых. Три вида включены в список Конвенции по международной торговле видами дикой фауны и флоры, находящимися под угрозой уничтожения (CITES). Однако остановить незаконный промысел, несмотря на декларируемые усилия, не удаётся, т. к. браконьерство является частью системной социальной проблемы. За счёт гибкой системы регулирования удалось предотвратить деградацию запасов горбуши тихоокеанского побережья Канады, где установился баланс «промысел - пополнение». Падение уловов привело к уменьшению (из-за низкой рентабельности промысла) количества купивших лицензию на лов рыболовецких судов. Заключение. Данные о динамике в системе «запас - пополнение - промысел» для популяций рыб почти всегда неполны, приблизительны и труднодоступны для стороннего исследователя. В отечественной практике редко можно говорить о наличии достоверных данных об изменении одновременно и запаса, и промыслового усилия. В ситуации с осетровыми Каспийского моря данные о промысловом усилии неизвестны. Сведения о современных размерах нелегального промысла крайне противоречивы. В работе [14] были опубликованы данные о динамике вылова, количестве орудий лова и состоянии запасов основных промысловых популяций Горьковского водохранилища. Сведения подтверждают динамику развития «оптимистичного» сценария. После начавшегося перелова численность популяции резко снижается и далее «псевдостабилизируется». В случае снижения промыслового усилия происходит быстрый рост численности видов рыб, т. к. эффективность воспроизводства из-за снижения конкуренции на ограниченных пригодных нерестилищах в водохранилищах достаточна. Промысел может использовать ∗ столь продуктивное состояние популяции, однако оно близко к R1 , и из-за кажущего незначительным превышения вылова может вновь произойти истощение ресурсов популяции на длительный период. Полученные с применением математических моделей выводы практически значимы при согласовании в наблюдениях или экспериментах. Для математической биологии подтверждение теории всегда было более проблематичным, чем для физики. Практикам управления эксплуатацией требуются средства, построенные именно в соответствии с особенностями жизненного цикла их промыслового объекта. В дальнейших исследованиях планируется описать в непрерывно-событийной модельной структуре влияние стадийности развития на эффективность рыбоводных мероприятий и рассмотреть сценарии модификации технологией искусственного воспроизводства биоресурсов. Предлагаемые методы вычислительного моделирования развиваются для системного обобщения накапливаемых экологами сведений о резких изменениях в биосфере.

About the authors

Andrey Yu Perevaryukha

St. Petersburg Institute for Informatics and Automation of RAS

Email: madelf@pisem.net
(Cand. Technical (Computer) Sci.), Senior Researcher, Lab. of Applied Informatics 39, 14-ya liniya V. O., St. Petersburg, 199178, Russian Federation

References

Supplementary files