Баканев С. В. Моделирование популяционной динамики камчатского краба на основе байесовского подхода // Принципы экологии. 2012. № 3. С. 4–23. DOI: 10.15393/j1.art.2012.1441

Выпуск № 3

Аналитический обзор

pdf-версия статьи

595.384.2-152.6.08(268.45)

Моделирование популяционной динамики камчатского краба на основе байесовского подхода

Введение

Одним из эффективных методов выявления механизмов функционирования популяции в ходе ее адаптации к новым условиям обитания является моделирование процессов пополнения, роста и убыли (Drury et al., 2007). Закономерности формирования популяционной структуры вида могут быть выражены количественно через параметры взаимодействия этих популяционных процессов. Несмотря на интенсивную двадцатилетнюю историю исследований камчатского краба Баренцева моря, многие аспекты оценки его популяционных параметров, в том числе и важные для промысла, остаются малоизученными. Моделирование индивидуального и группового роста ракообразных, в отличие, например, от рыб, сопряжено с трудностью определения точного возраста. Камчатский краб растет в течение короткого периода после линьки, т. е. индивидуальный рост особей имеет сложный прерывистый характер. Частота линьки зависит от пола, возраста и условий обитания краба. Традиционные модели с использованием численности поколений становятся трудноприменимыми (Smith, Addison, 2003). Кроме того, камчатский краб для Баренцева моря – инвазийный вид, интродукция которого была осуществлена в 60-е годы прошлого столетия. При взрывном характере численности, который часто наблюдается у акклиматизированных животных, моделирование динамики численности и оценка экологической емкости среды для интродуцента имеет существенные неопределенности и затрудняет принятие допущения о динамическом равновесии популяции, необходимом в традиционном моделировании системы «окружающая среда – запас – промысел». В период открытия коммерческого промысла камчатского краба (2004–2007 гг.) одной из главных проблем в определении продуктивности баренцевоморской популяции камчатского краба являлся его незаконный промысел, который не поддавался точной оценке. По экспертным оценкам российских исследователей, незаконный экспорт, нелегальный вылов и браконьерство, наряду с несовершенной системой контроля промысла, не позволяли оценить реальные объемы добычи краба и, соответственно, затрудняли определение уровня эксплуатации и продуктивности популяции. В процессе разработки модельных подходов к динамике численности камчатского краба в Баренцевом море был выработан алгоритм аналитической оценки, который учитывает все вышеизложенные аспекты, однако сами подходы по сей день подвергаются изменениям. Причины таких изменений связаны, прежде всего, не со стремлением усовершенствовать математический аппарат процедуры, а с необходимостью учитывать меняющееся качество и объем исходных данных, получаемых в ходе ежегодных исследований. Сокращение объемов финансирования исследований и, как следствие, изменения в методике проведения исследовательских съемок стали причиной пересмотра существующих процедур оценок в последние годы. В условиях плохой информационной обеспеченности исследований величины параметров зачастую находят, используя алгоритмы на основе формулы Байеса, когда в качестве исходной информации берутся не только данные наблюдений, но и априорное (предварительное) знание о параметрах модели (Bayes, 1763). История применения подходов на основе теоремы Байеса в популяционной биологии насчитывает порядка 20 лет (Burgner, 1991; Angel et al., 1994; Schnute, 1994; Punt, Butterworth, 1996; Kinas, 1996; Meyer, Millar, 1999; Patterson, 1999). В настоящее время такие подходы широко используются для описания состояния и прогнозирования динамики запасов морских млекопитающих, рыб и беспозвоночных. Модели, использующие теорему Байеса, были успешно адаптированы для оценки запасов северной креветки, синего краба, канадского лосося, южноафриканского анчоуса, усатых китов, тихоокеанских тунцов, палтуса и сельди (Butterworth, Punt, 1995; Schnute et al., 1998; Meyer, Millar, 1999; Millar, Meyer, 2000). Оценки состояния запасов и прогнозы, выполненные при помощи таких моделей, используются в различных международных организациях, участвующих в формировании рекомендаций по эксплуатации гидробионтов: ФАО, АНТКОМ, НАФО, ИКЕС, НЕАФК (Annala, 1993;  Adkison and Peterman, 1996; McAllister and Kirkwood, 1998; Punt, Hilborn, 2001; Hvingel and Kingsley, 2006;  ICES, 2008). В настоящей работе описан опыт моделирования динамики популяции камчатского краба в Баренцевом море на основе байесовского подхода. Дается краткий обзор методов когортного и продукционного анализа популяционных параметров с учетом экологических особенностей этого ракообразного и методов сбора информации. Рассматривается методика применения байесовского подхода на примере оценки параметров продукционной модели. 

Аналитический обзор

 Модели популяционной динамики

Для оценки биологических показателей популяции камчатского краба в Баренцевом море в 2003 г. была осуществлена первая попытка моделирования  динамики его численности. Модель описывала элементарные биологические процессы и базировалась на традиционных биостатистических методах, разработанных для оценки эксплуатируемых популяций рыб. Рассматривая динамику численности одного поколения, профессор Ф. И. Баранов вывел в 1918 г. соотношение, на которое опирается в настоящее время большинство биостатистических моделей (Баранов, 1918):

, [1]

где N_t – численность поколения в начальный момент времени t, N_t+1 – численность в момент времени t + 1, e – основание натуральных логарифмов. Показатель степени Z называется коэффициентом мгновенной общей смертности (Рикер, 1979; Шибаев, 2007). Смысл коэффициента: за элементарный (очень маленький) промежуток времени dt численность животных уменьшается на величину dN, равную Z-й части от фактической численности N или dN / dt = –ZN. Для эксплуатируемых популяций гидробионтов Z определяется изъятием животных промыслом (коэффициентом мгновенной промысловой смертности, F) и всеми остальными причинами, кроме промысла (то есть коэффициентом мгновенной естественной смертности, M) :  Z = F + M.

Спустя 40 лет J. A. Gulland (1964) и G. I. Murphy (1965) предложили метод анализа структуры уловов эксплуатируемых популяций рыб на основе подхода Ф. И. Баранова. Использовались два уравнения. Первое отражает взаимосвязь численностей смежных возрастных групп:

, [2]

Второе уравнение позволяет рассчитать величину улова в поштучном выражении, получаемого от каждой возрастной группы:

, [3]

где N_{t, a} , N_{t+1, a+1} – начальная численность возрастной группы a в году t и в следующем t + 1 году, когда возраст рыбы увеличится на единицу; C_F – улов возрастной группы a в году t; F_t,a, M_t,a – коэффициенты мгновенной промысловой и естественной смертности возрастной группы a в году t. Они могут быть зависимыми или независимыми от возраста, а величина F может изменяться по годам в связи с изменениями интенсивности промысла. Уравнения служат базовыми для группы методов, широко используемых в настоящее время в промысловой ихтиологии, которая называется виртуально-популяционным анализом ( Murphy, 1965; Pope, 1972; Jones, 1974; Рикер, 1979).

При виртуально-популяционном анализе, как правило, начальная численность N_t,1 определяется величиной пополнения R_t,1, то есть количеством особей, пополнивших в году t промысловый запас (Рикер, 1979). Обозначим количество животных, погибших по естественным причинам, символом C_t,a, тогда схему эксплуатируемой популяции рыб, согласно уравнениям [2] и [3], можно представить графически (рис. 1).

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

Fig.1. Block diagram of modeling the  dynamics of the size of an exploited population (solid arrows – the transition of survivors into the next age group, the dashed arrows – the loss of animals from fishing and natural mortality)

 

Виртуально-популяционный анализ динамики численности промысловых ракообразных на основе возрастных когорт не получил широкого применения в силу сложности определения точного возраста животного (ICES, 2001a; ICES, 2001б). Однако с начала 90-х годов прошлого столетия, основные принципы изменения численности поколений, разработанные Ф. И. Барановым, стали с успехом применяться в моделировании динамики численности крабов и креветок (Collie, 1991; Kruse and Collie, 1991; Zheng et al., 1995; Quinn et al., 1998; Cadrin, 2000). Модели этого типа базировались на изменениях численности не поколений, а размерных групп (LBA, length-based analysis) и учитывали линьку и стохастический рост животных.

Для популяции камчатского краба Баренцева моря модель LBA была реализована в 2003 году и стала использоваться с 2006 года, когда временной ряд наблюдений составил 13 лет (Баканев, 2006). Рассмотрим принцип моделирования динамики численности эксплуатируемой популяции камчатского краба на основе модели LBA, описанной в работе J. Zheng и др. (1995). Численность краба моделировалась раздельно по размерным группам l = 1, 2, 3, … i и по полу. Каждая размерная группа l в году t состояла из особей с новым (N_l,t) и старым (O_l,t) карапаксом. В течение года t крабы как с новым, так и старым панцирем могут линять с вероятностью m_l,t или с вероятностью 1 – m_l,t, не линяя, оставаться в этой же размерной группе.

Вероятность линьки для краба m_l,t была представлена логистической функцией (Balsiger, 1974):

 , [4] 

где a_t и b_t – параметры, ^l – средняя длина размерной группы l.

Нелиняющие крабы N_l,t и O_l,t в году t + 1 становились O_l,t+1 в таком соотношении:

 , [5] 

Линяющие крабы в t году – m_l,t( N_l,t + O_l,t), в году t +1 в зависимости от длины прироста попадают в последующие размерные группы и становятся крабами с новым карапаксом, принадлежащими размерным группам 1, 2, 3, …, L. Вероятность перехода в определенную группу зависит от параметров прироста краба за одну линьку. Рост на индивидуальном уровне характеризуется длиной прироста особи. Длина прироста особи на групповом уровне варьирует как внутри размерной группы, так и между группами. Средний прирост особи камчатского краба для размерной группы l (G_l,) обычно выражали линейной функцией от средней длины (^l) размерной группы l перед линькой краба (Weber, Miahara, 1962):

, [6]

где a и b – параметры линейной функции.

Отклонения прироста особей размерной группы l за линьку от уравнения регрессии описывали через вероятностное распределение. Плотность вероятности величины прироста выражали через гамма-распределение, форма которого, в зависимости от значений параметров, охватывает все возможное многообразие изменений в распределении прироста (Sullivan et al., 1990):

, [7]

где x – величина прироста во время линьки, a_l и b – параметры. Средняя x задавалась произведением a_l · b и равна G. Так a_l = G / b, тогда прирост выражался посредством двух параметров G и b. Ожидаемая пропорция крабов размерной группы l, переходящих после линьки в последующую группу l', равна интегралу функции гамма-распределения на интервале группы l':

 ^, [8]

где ^l – средняя длина группы l, l₁ – начало интервала группы l', l₂ – конец интервала группы l'. Для последней размерной группы L, P_LL = 1.

Пополнение (R) моделируемой популяции камчатского краба американские исследователи описали с помощью двух параметров: 1) количество молодых крабов, пополняющих популяцию в году t, R_t и 2) параметр U_l, определяющий стохастическое изменение процесса пополнения (Sullivan et al., 1990; Zheng et al., 1995). Пополнение каждой размерной группы R_l,t выражалось через уравнение:

, [9]

где U_l описывалось гамма-распределением, таким же как в уравнениях [7] и [8] с параметрами a_r и b_r и. Средний размер пополнения задавался произведением a_r · b_r, поэтому требовалось оценить только параметр b_r, при условии, что был известен средний размер пополнения.

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

, [10]

где N_t – вектор L размерных групп, показывающий общее количество животных с новым карапаксом в году t. Расчет N_t+1 обеспечивает матрица переходов (transition matrix), содержащая параметры стохастического роста (m, P) и пополнения (R):

, [11]

где A_l,t – временная переменная, показывающая  остаточную численность краба с новым и старым карапаксом в размерной группе l с учетом вылова, и равная:

, [12]

где N_l,t и O_l,t – численности крабов с новым и старым карапаксом в размерной группе l в году t соответственно; C_l,t – вылов крабов (экз.) размерной группы l в году t; M_l,t – коэффициент мгновенной естественной смертности для крабов размерной группы l в году t и y_t – временнáя задержка от момента съемки до средней даты периода промысла. Расчет переменной A_l,t основывается на подходе Ф. И. Баранова, описанном уравнениями [2] и [3]. Он также схематично представлен на рисунке 1.

Основная сложность реализации методов, основанных на когортном анализе, заключается в том, что для исходных систем уравнений трудно получить однозначное решение, особенно в случае низкого качества входных данных (Бабаян, 2000). Так, для баренцевоморской популяции камчатского краба короткий ряд наблюдений, относительно низкие уловы крабов в ходе ежегодных съемок, недооценки численности вследствие активных миграций акклиматизанта, неопределенности в оценке общего вылова, а также методические ошибки в ходе становления нового ряда наблюдений вносили существенные помехи в оценку параметров запаса и затрудняли выявление закономерностей его динамики.

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

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

На основе модели LBA в 90-е годы прошлого столетия был разработан целый ряд моделей, успешно использующих возможности анализа размерных рядов (Collie, 1991; Kruse and Collie, 1991; Zheng et al., 1995; Quinn et al., 1998; Cadrin, 2000). Одна из них, CSA (catch survey analysis), является частным случаем LBA, когда численность популяции мала для получения качественных данных по размерному составу. Модель была разработана Зенгом и др. (Zheng et al., 1997) для небольших запасов краба с относительно низкими уловами в ходе ежегодных съемок. Обе модели принципиально сходны. Отличие состоит лишь в том, что в модели CSA используются более укрупненные размерные группы, а в расчетах приемлемо использовать данные по размерному составу лишь половозрелых самцов. В модели CSA деление на размерные группы, как правило, основывается на особенностях биологии и промысла дальневосточных крабов семейства Lithodidae, в которое входит и камчатский краб. Промысловая часть популяции в году t, состоящая из крупных самцов, делится на две группы: рекруты (RE_t) и пострекруты (PO_t). RE_t – это численность крабов с новым карапаксом, пополнивших промысловый запас в году t. Рекруты принадлежат размерной группе, начало которой равно минимальному промысловому размеру, а длина интервала размерной группы соответствует средней длине прироста рекрута за одну линьку. PO_t– численность крабов, пополнивших промысловый запас в году t – 1 и ранее. Пострекруты это все промысловые крабы со старым карапаксом, а также особи с новым панцирем, имеющие размеры большие, чем рекруты. Пополнением промыслового запаса считаются пререкруты (PR_t) – особи, которые в следующем году станут рекрутами. Пререкруты принадлежат размерной группе, конец интервала которой граничит с минимальным промысловым размером, а длина интервала соответствует средней длине прироста пререкрута за одну линьку. Динамику моделируемой популяции можно представить в виде уравнений:

^, [13]

^, [14]

где m_PR – вероятность линьки для пререкрутов; G_PR,RE, G_PR,RO – параметры роста, соответствующие доле крабов, которые переходят из группы PR_t в группу RE_t+1 и из группы PR_t в группу PO_t+1; M – коэффициент естественной смертности; C_t – вылов крабов в году t; y_t – временнáя задержка от момента съемки до средней даты периода промысла.

Адаптация опыта американских ученых и результаты наших исследований показали хорошие перспективы для использования стохастической версии модели CSA как инструмента для оценки динамики запаса, отслеживания реакции популяции на промысел, оценки ОДУ и ориентиров управления для баренцевоморского запаса камчатского краба (Баканев, 2003, 2008). Результаты оценок популяционной динамики камчатского краба в Баренцевом море по модели CSA показали лучшую согласованность исходных индексов численности и модельных значений по сравнению с результатами использования модели LBA. Кроме того, относительная простота и меньшее количество параметров модели CSA делают ее менее чувствительной к качеству исходной информации.

Описанные ранее модели LBA и CSA, основанные на динамике численности размерно-возрастных групп, требуют как минимум наличия исходных данных о размерно-возрастной структуре популяции на протяжении достаточно продолжительного временного ряда. Такая исходная информация об эксплуатируемой популяции, как правило, собирается в результате ежегодного мониторинга промысла и выполнения научно-исследовательских съемок. Нарушение методик выполнения съемок в последние годы привело к поиску альтернативных методов оценок динамики численности популяции камчатского краба. При дефиците биологической и промысловой информации моделируемая численность может быть выражена упрощенно через изменение производительности промысла и величине вылова (Бабаян, 2000; Баканев, 2009).

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

Для баренцевоморской популяции камчатского краба продукционная модель была адаптирована в 2006 году с открытием полномасштабного коммерческого промысла краба. Дискретная форма записи модели при отсутствии промысла имеет вид (Schaefer, 1957):

^, [15]

где B_t – биомасса животных в году t, K – максимальная возможная биомасса при отсутствии промысла (емкость среды), r – коэффициент мгновенного популяционного роста (внутренняя скорость роста) в отсутствие плотностной регуляции. Согласно свойству этого логистического уравнения максимальные продукционные способности запаса наблюдаются при биомассе, равной K / 2.

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

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

Другой альтернативой размерно-возрастным моделям в настоящее время становится модель истощения Лесли (Leslie depletion model; Haddon, 2001). Модель основана на динамике снижения производительности промысла в течение промыслового сезона и используется для оценки численности на начало промыслового сезона с учетом темпов снижения.

Учитывая относительно короткий период промысла (2–3 месяца), делается допущение, что величина естественной смертности за период промысла незначительна, и ей можно пренебречь. Тогда уравнение динамики численности на акватории промысла будет выглядеть так: N_t+1 = N_t + C_t. Численность на начало промыслового сезона (N_t+1) при этом будет равна сумме численности на конец сезона (N_t) и кумулятивному вылову (C_t) между временем 1 и t. При этом N может быть выражен через улов на усилие или производительность (U) и коэффициент улавливаемости: U = qN, где U – улов на усилие, q – коэффициент улавливаемости. С учетом коэффициента улавливаемости, уравнение динамики численности на акватории промысла будет выглядеть так: U_t+1 = U_t + qC_t. Уравнение имеет вид линейной регрессии, с помощью которого может быть вычислен q и, соответственно, N. Для оценки численности камчатского краба в начале промыслового сезона в Баренцевом море используются данные промысловой статистики 2007–2012 гг., с учетом того, что в этот период промысел был локализирован в одних и тех же районах, а степень неучтенного вылова значительно снизилась. Характер снижения производительности с учетом возрастающего вылова к концу промыслового сезона используется для расчета численности запаса на акватории промысла в период этого промыслового сезона.

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

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

1) построение модели исследуемого явления;

2) оценка основных характеристик (параметров) модели по данным наблюдения.

Основы первого этапа, описанные в предыдущей главе, заложены к середине прошлого века и не претерпели существенных изменений (Рикер, 1979, Одум, 1986). Второй этап опирается на методы математической статистики, которые в настоящее время продолжают активно развиваться (Punt and Hilborn, 1997; Schnute et al., 1998).

Стохастическая природа оцениваемых параметров может быть описана с использованием различных алгоритмов идентификации, наиболее распространенным из которых является метод максимального правдоподобия (Hilborn and Walters, 1992). Точность получаемых оценок параметров зависит от полноты и качества исходной информации. Фрагментарность данных, ошибки измерений, шум природных процессов в совокупности с короткими рядами наблюдений могут значительно исказить результаты оценок и, соответственно, понимание реальных биологических явлений (Васильев, 2001).

В условиях плохой информационной обеспеченности оценки параметров зачастую находят, используя алгоритмы на основе формулы Байеса, когда в качестве исходной информации берутся не только данные наблюдений, но и априорное (предварительное) знание о параметрах модели (Bayes, 1763). Подход основан на попытке начать статистический вывод с некоторых исходных предположений (догадок) о вероятностном распределении неизвестных параметров. Значения модельных параметров можно установить априорно на основе оценок, полученных, например, для одних и тех же видов, но из разных районов, или для схожих видов из одного района. Используя теорему Байеса, окончательные (апостериорные) значения параметров оцениваются с учетом как данных наблюдений, так и предварительно заданных значений параметров:

^, [16]

Здесь распределение P(H) называется априорным распределением вероятностей возможных значений параметра или параметров H (это распределение принимается прежде, чем получены статистические данные); data – статистические данные, полученные в ходе эксперимента (наблюдений, опытов), которые используются в модели; P(H | data) – условная апостериорная вероятность, количественно оценивает насколько модель соответствует гипотезе H в условиях, характеризуемых данными наблюдений (data); P(data | H) – функция правдоподобия, используемая в традиционном статистическом оценивании параметров, которая определяет вероятность того, что полученные данные в ходе эксперимента соответствуют гипотезе H .

В рамках байесовского подхода, гипотезы (параметры) трактуются как случайные величины и описываются с использованием вектора непрерывных значений параметров вместо набора дискретных значений. Вероятности в этом случае будут выражаться через плотности вероятностей, а знаменатель P(data) в уравнении [16] согласно формуле полной вероятности будет равен

.

Функцию Байеса можно представить иначе, в виде трех составляющих:

^, [17]

где 1/P(data) некая постоянная величина по отношению к гипотезе. Поэтому апостериорное распределение пропорционально произведению функции правдоподобия на априорное распределение или:

^, [18]

где символ µ  означает пропорциональность.

 В терминах статистики байесовский подход базируется на следующих четырех компонентах. Априорная (предварительная) информация (PRIORS) и данные (DATA), полученные в ходе эксперимента, объединяются с помощью некоего метода (method) для получения апостериорного (POSTERIORS) знания. Эти компоненты могут быть представлены в виде:

^, [19]

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

На первом этапе первой итерации происходит расчет численности животных (data) по биологической модели (BIO) с использованием стартовых (начальных) значений параметров и переменных (INITS). Биологическая модель в нашем случае представляет собой одну из четырех, описанных ранее: LBA, CSA, продукционную или модель истощения. Стартовые значения генерируются встроенным модулем программы, а затем, если необходимо, корректируются с учетом биологически правдоподобных значений. На первом этапе второй и последующих итераций берутся величины параметров и переменных, рассчитанные в ходе предыдущей итерации.

На втором этапе модельные и эмпирические (DATA) значения численности включаются в функцию правдоподобия (LIKELIHOOD) и вычисляется (like) мера их согласованности, то есть определяется вероятность того, что эмпирические значения могут быть получены с использованием данной биологической модели с заданными стартовыми значениями параметров и переменных.

На третьем этапе по формуле Байеса (BAYES, уравнение [18]) количественно оценивается способность модели с заданными параметрами (PRIORS) генерировать эмпирические значения DATA. В процессе расчетов происходит настройка параметров (PRIORS) в условиях, характеризуемых данными наблюдения (DATA), и рассчитываются новые апостериорные значения модельных параметров и переменных. Оцененные показатели параметров и рассчитанные модельные переменные, характеризующие модель в данный момент, играют роль начальных условий для следующей итерации.

Данная процедура расчета является заключительным этапом одного итерационного шага или звена в цепи Маркова (MCMC) и реализована с использованием метода Гиббса. В процессе итераций генерируются (stat) выборки распределений значений параметров и расчетных переменных. Статистический анализ таких выборок, в который входит расчет средних, ошибок средних, стандартных отклонений, медиан и доверительных интервалов параметров, является результатом прогона модели. Окончание итерационного процесса (прогона) происходит, когда ошибка среднего искомого параметра, рассчитанная в цепи Маркова, составит менее 5 % его стандартного отклонения (Spiegelhalter et al., 2000). 

 

Рис. 2. Принципиальная блок-схема оценки параметров биологических моделей (показаны два итерационных шага). DATA – эмпирические данные; INITS – стартовые значения параметров и переменных; PRIORS – априорные распределения параметров; POSTERIORS – апостериорные распределения параметров и переменных; BIO – биологическая модель (уравнение популяционной динамики); LIKELIHOOD – функция правдоподобия; BAYES – формула Байеса; input – ввод данных; data – расчет численности животных по биологической модели; like – расчет правдоподобия; MCMC – переход к следующей итерации; stat – расчет апостериорных распределений параметров и переменных

Fig. 2. Schematic block diagram of the parameter estimates of biological models (shown in two iterative steps). DATA – empirical data; INITS – starting values ​​of parameters and variables; PRIORS – a prior distribution of the parameters; POSTERIORS – the posterior distribution of parameters and variables; BIO – biological model (thq equation of population dynamics); LIKELIHOOD – likelihood function; BAYES – Bayes' formula; input – data input; data – calculation of the number of animals according to a biological model; like – the likelihood calculation; MCMC – the transition to the next iteration; stat – the calculation of posterior distributions of the parameters and variables

 

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

Во-первых, рассчитывались остатки, т. е. разности между логарифмами фактических индексов численности по съемке и рассчитанных индексов по модели. Чем лучше модель согласуется с данными, тем меньше величина остатков. Суммарная статистика отклонений (RSS, от residual sum of squares) показала степень смещения модельных значений от наблюдаемых всего временного ряда:

^,

где N_t – наблюдаемый индекс численности по съемке, N'_t – расчетный индекс численности по модели.

Во-вторых, адекватность модели анализировалась путем сравнения распределения каждой фактической величины со своим апостериорным прогнозным распределением (Gelman and Rubin, 1992). При этом после настройки параметров в модели имитируется набор наблюдаемых данных. Вероятности имитируемых значений фактических величин data_rep имеют распределение:

^, [20]

В выражении [20] показатель P(data_rep |q ) представляет собой выборку наблюдений из распределений, рассчитанных в модели – P(q | data_obs).

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

^, [21]

где I(x) = 1, если x истинно, и I(x) = 0, если x ложно. Величины, близкие к 1 или 0 в векторе p-значений, показывают, что набор фактических значений в значительной мере отклоняется от медианы апостериорного распределения. Величина p.value близкая к 0,5 показывает, что расчетная величина data_j^rep приблизительно в 50 % случаев больше входного значения и в 50 % – меньше, то есть отклонения не имеют систематического смещения.

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

Реализация модельных алгоритмов и статистическая обработка полученных результатов выполняется в программном продукте OpenBUGS, разработанном в Медицинском исследовательском центре в Кембридже (Spiegelhalter et al., 2000; OpenBUGS, 2009). Программа снабжена руководством пользователя, в котором подробно описан байесовский алгоритм пошаговой реализации математической модели, оценки параметров и диагностики.

 

Реализация байесовского подхода в оценке параметров модели

 

Рассмотрим методику применения байесовского подхода на примере оценки параметров продукционной модели Шефера, описанной выше  уравнением [15]. С добавлением показателя вылова модель приобретает вид (Schaefer, 1954):

,[22]

где B_t – численность в году t, C_t – вылов в году t, K – максимально возможная численность популяции и r – коэффициент мгновенного популяционного роста. Ненаблюдаемая переменная B_t может быть выражена через наблюдаемый показатель относительной численности (I_t), рассчитанный, например, по результатам исследовательской съемки:

^, [23]

где q – коэффициент пропорциональности, в данном случае коэффициент улавливаемости, a  – остаточная погрешность, имеющая логнормальное распределение (Haddon, 2001).

В системе уравнений [22] и [23] нужно оценить четыре параметра: r, K, q, , а также начальное значение численности B₀. Традиционно в популяционной биологии оценка параметров осуществляется либо методом наименьших квадратов, либо методом максимального правдоподобия (Hilborn and Walters, 1992). Как показано выше, функция правдоподобия  также входит в состав формулы Байеса [17] и выглядит так:

, [24]

Функция правдоподобия определяет вероятность получения данных наблюдений (data) при определенном наборе параметров, значения которых соответствуют гипотезе H, то есть r, K, q, и B₀. Метод максимального правдоподобия заключается в подборе таких значений параметров, при которых вероятность  P(data|r, K, q,, B₀) стала бы максимальной. Чем выше вероятность, тем лучше модель описывает данные наблюдений, а следовательно, динамику численности популяции. При дефиците исходной информации результаты настройки модели могут привести к неверным выводам о численности животных. Например, если дать нереально высокую величину K или q и тем самым значительно исказить расчетную численность. Для предотвращения таких артефактов необходимо внести дополнительную информацию, то есть некое априорное (доопытное) знание.

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

Отображая соотношение реальной величины численности и наблюдаемого индекса относительной численности, он может быть предварительно оценен посредством полевых экспериментов. Являясь коэффициентом пропорциональности для индекса численности, оцененной в ходе исследовательской съемки, этот показатель определяется, главным образом, уловистостью орудий лова. Доля пойманных животных из числа находящихся в зоне отлова, естественно, может принимать значения от 0 (никто не пойман) до 1 (отловлены все особи). На основании этой информации эксперт может допустить, что действительное значение коэффициента q равновероятно может находится в диапазоне от 0 до 1.

Допустим, что у эксперта имеется некая дополнительная информация в виде литературных данных о коэффициенте уловистости аналогичного трала по отношению к животным исследуемого нами вида. В результате он может прийти к выводу о том, что интервальная оценка коэффициента уловистости трала должна лежать, например, в пределах 0.4–0.6. Кроме того, из многолетнего опыта эксперта и литературных данных известно, что коэффициент уловистости варьирует от траления к тралению и зависит от множества факторов (параметров раскрытия трала, типа грунта, размеров животных и пр.). На основании этой информации, а также используя соображения о характере варьирования этого показателя, исследователь имеет некое мысленное представление о возможной величине этого параметра до начала моделирования процесса. Ожидается, что при использовании традиционного метода максимального правдоподобия и данных наблюдений, величина q окажется, скорее всего, близкой к 0.4–0.6. В тоже время маловероятно, что значение этого коэффициента будет близким, например, к 0 или к 1.

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

Построение априорного распределения целесообразно начинать с выбора типа распределения. Однако на практике не всегда имеется исчерпывающая информация о том, какому закону подчиняется вероятностное распределение того или иного параметра. В настоящее время проблема выбора априорных распределений для параметров моделируемых биологических процессов освещена подробно (Kass, Wasserman, 1996; Punt and Hilborn, 1997). В данном случае, учитывая литературные данные, предполагаем, что распределение подчиняется нормальному закону с модой в интервале 0,4–0,6 и границами от 0 до 1. Графически такое распределение изображено на рисунке 7, а математически принято записывать как:

q~dnorm(0,50;50), [25]

где dnorm – нормальное распределение, а символ "~” обозначает “распределяется как”.

Параметры нормального распределения mu (среднее) и tau (точность, показатель обратно пропорциональный дисперсии, то есть ) указаны в скобках. При включении априорного распределения q в формулу Байеса [17] рассчитывается апостериорное распределение этого параметра. Иными словами наша предварительная, или априорная, оценка параметра q вновь оценивается с учетом данных наблюдений, априорное распределение корректируется и строится окончательное апостериорное распределение (рис. 3). На этом процесс оптимизации или настройки параметра q заканчивается. В данном случае была получена уточненная оценка q с медианой, равной 0,38.

 

Рис. 3. Распределение вероятности параметра q (черная линия – априорное, синяя – апостериорное)

Fig. 3. Probability distribution of the parameter q (black line – a priori, blue – a posteriori)

 

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

 

 

Рис. 4. Распределение вероятности параметра q (черная линия – априорное, синяя – апостериорное) при информативных (А) и малоинформативных (Б) данных наблюдений

Fig. 4. Probability distribution of the parameter q (black line – a priori, blue – a posteriori) for informative (A) and uninformative (B) observations+

 

Наиболее важным и одновременно сложным является вопрос выбора априорного распределения параметров. Для определенных параметров выбор может основываться в некоторых случаях на субъективном мнении эксперта. В идеале это мнение должно быть основано на фактах. Эти факты могут включать все те сведения, которые обычно не входят в формальный анализ количественных данных, а используются как информация, полученная вне эксперимента. Такая «внешняя» информация может включать:

1) результаты исследований по той же тематике;

2) результаты исследований, в которых изучались сходные биологические механизмы;

3) результаты лабораторных исследований, посвященные природе изучаемого явления;

4) результаты исследований, если эти результаты могут иметь ту же природу;

5) сведения о промежуточных исходах, которые наблюдались в данном эксперименте и свидетельствуют в пользу предложенной гипотезы;

6) информацию, полученную у других популяций со схожими биологическими и промысловыми параметрами.

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

Согласно формулам [20] и [21] для настройки всего моделируемого процесса должны быть заданы априорные распределения не только для параметра q, но и для r, K, и B₁. Далее на основе априорных распределений пяти параметров моделируется общее априорное распределение. Затем, используя формулу Байеса и данные наблюдений, рассчитывается целевое или общее апостериорное распределение, статистические показатели которого являются искомыми оценками параметров модели.

Рассмотренный пример байесовского подхода на примере продукционной модели представляет собой принципиальную схему оценки параметров модели. Реализация формулы Байеса, когда все параметры и расчетные переменные трактуются как случайные величины, основана на поиске решения системы интегральных уравнений. Процедура поиска решения связана с большим количеством вычислений и базируется на итеративном методе с применением генератора случайных чисел. В настоящей работе используется подход, который в англоязычной литературе известен как Markov Chain Monte Carlo (MCMC; Congdon, 2001). Этот алгоритм вычисления основан на построении цепи Маркова. Процесс, протекающий в системе, называется марковским, если для каждого момента времени которой вероятность любого состояния системы в будущем определяется только состоянием системы в настоящий момент и не зависит от того, каким образом система пришла в это состояние (Вентцель, 2003). Реализация цепи Маркова в численном виде происходит с помощью модели Гиббса (Gibbs sampler). Она обеспечивает способ выборки из совместных распределений многомерных переменных с помощью применения такого положения: для получения выборки из совместного распределения делаются многократно выборки из его представленных одномерных условных распределений. Процесс реализации MCMC с помощью модели Гиббса подробно описан в литературе (Gilks et al., 1996; Meyer, Millar, 2000).

Заключение

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

Оценка состояния запаса камчатского краба в Баренцевом море основывается на реализации традиционных методов определения продуктивности популяции. В основе оценки лежат основные постулаты теории рыболовства, в соответствии с которыми выполняется расчет популяционных параметров. В то же время биологические и экологические особенности камчатского краба и история его появления в Баренцевом море вынуждают отходить от алгоритмов, традиционно используемых при оценке запасов гидробионтов. Во-первых, трудность в определении возраста ракообразных не позволяет использовать преимущества моделей, основанных на оценке численности поколений. Во-вторых, характер роста численности искусственно созданной популяции затрудняет использование принципов управления, основанных на динамическом равновесии системы «окружающая среда – запас – промысел». В-третьих, на этапе начала промысловой эксплуатации запаса возникают трудности интерпретации связи биологических параметров популяции с промысловыми показателями. Кроме того, прямые наблюдения по результатам съемок и промысла несут значительное количество шума и зачастую весьма неопределенны, что связано как с методическими проблемами, так и с экологическими особенностями вида.

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

Библиография