Регулирование многовидового рыболовства на основе математического моделирования (стр. 39 )

M=a + Can·P ,

где Can - постоянная скорость смертности от каннибализма, a – коэффициент естественной смертности от других причин.

Чаще всего вместо численности в эту формулу подставляют биомассу нерестового запаса (SSB). Е. Нильсен с соавторами (Nilssen et al., 1994) выбирают в качестве индекса численности хищников P северо-восточной арктической трески суммарную численность возрастных групп 3 и 4 трески. И то и другое не совсем точно, поскольку и половозрелые особи всех возрастов и неполовозрелые особи возрастных групп 3 и 4 трески потребляют свою молодь.

Уравнение модели Рикера можно представить в следующем виде: численность поколения в возрасте tr (в нашем случае, в возрасте трех лет), обозначаемая через R и называемая пополнением, равна

. (7-1)

Здесь N(0) – начальная численность поколения, например, количество отложенной икры. Если константа tr включена в параметры a и Can, а N(0)=mSSB(Y) (где m=const, а SSB(Y) – биомасса нерестового запаса в год рождения поколения Y), то модель имеет вид

(7-2)

Рассмотрим несколько предположений или упрощений, сделанные У. Рикером при разработке модели запас-пополнение:

·  Каннибализм действует одинаково на всех стадиях формирования пополнения, от возраста 0 до tr , то есть параметр Can=const.

·  Начальная численность поколения N(0), а именно общее количество отложенной икры или популяционная плодовитость, пропорциональна биомассе нерестового запаса.

·  Численность хищников (каннибалов) P пропорциональна биомассе нерестового запаса SSB .

Параметры уравнения модели Рикера (7-1) оцениваем регрессионным методом, минимизируя сумму квадратов отклонений двух временных рядов – численности пополнения в возрасте трех лет, полученной на Рабочей группе ИКЕС по арктическому рыболовству AFWG (ICES, 2008) и оценок пополнения, полученных по модели (7-1). Первый ряд рассматриваем при этом как истинные значения пополнения, а второй – теоретические или модельные его значения. Выбраны поколения трески рождения с 1984 г. по 2000 г. Корреляционное отношение для этих двух рядов равно 0,33, что означает, что модель Рикера объясняет 33% изменчивости пополнения за этот период времени.

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

Второе предположение, сделанное при выводе уравнения Рикера, следует уточнить, поскольку количество отложенной икры зависит не только от величины нерестового запаса, но и от его возрастной структуры. Потому первый множитель уравнения (7-2) заменяем на величину популяционной плодовитости.

На успешность выживания пополнения по-видимому значительное влияние оказывает интенсивность и продолжительность притока атлантических вод, так как образующиеся при этом вихри вдоль основных потоков создают благоприятные условия для развития и кормового планктона и личинок трески (устное сообщение д. г.н. , ВНИРО). Чтобы формализовать это явление, введен показатель, рассчитываемый как количество месяцев в год рождения поколения, в которые наблюдались положительные температурные аномалии. Ясно, что этот показатель, обозначенный как N(Y), может меняться от нуля до 12. Делаем предположение, что этот индекс N(Y) характеризует установившийся приток атлантических вод через Кольский разрез, и с его увеличением коэффициент естественной смертности должен уменьшаться.

Уравнение модели теперь принимает вид:

, (7-3)

где Y – год рождения поколения,

Rmod (Y+3) –численность пополнения в возрасте трех лет (в 1000 экз.),

IndPF(Y) –индекс популяционной плодовитости, равный PF*10-12 (PF–population fecundity),

SSB(Y) – биомасса нерестового запаса в год Y, т.

Индекс притока атлантических вод N(Y) рассчитываем по температурным аномалиям (TA) на Кольском разрезе в слое 0-200 м в год рождения поколения.

Параметры a, Can и ny – константы, которые следует оценить.

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

На Арктической рабочей группе ИКЕС в 2008 г. (ICES, 2008) не было принято единого решения о величине несообщенного вылова трески в последние годы, потому там проводились две серии расчетов запаса с оценками несообщенного вылова, представленными норвежской стороной (этот расчет обозначен как IUU-Норв) и российской стороной (IUU-Росс). Потому параметры модели пополнения определены нами дважды – по результатам каждого из вариантов оценки запаса трески.

Входные данные для оценки параметров модели пополнения:

§ временные ряды R(Y+3) и SSB(Y), полученные с помощью когортной модели XSA с терминальным годом 2007 (два варианта прогонов модели с разными IUU);

§ значения PF (Y) рассчитаны по возрастному составу запаса, полученному в тех же двух прогонах, по методу М. Бондаренко и др. (2003), при этом взяты переменный по годам темп полового созревания (Yaragina, 2001) и возрастная зависимость индивидуальной плодовитости (Ярагина и др., 2003). Для всех лет доля самок в нерестовом запасе принята равной 0,5.

§ T(Y, m) – матрица среднемесячных температур в слое 0-200 м на станциях 3-7 Кольского разреза, полученная на сайте ПИНРО www. pinro. ru/labs/, по которой рассчитана матрица температурных аномалий относительно средних месячных температур за период с 1951г. по 2000 г., а затем и ряд N(Y).

Параметры модели пополнения оценены с помощью надстройки Solver в среде EXCEL для двух интервалов лет рождения поколений трески Y: 1984-2000 гг. и 1984-2003 гг.

В течение последних 10 лет отмечается сильно выраженный тренд в изменениях TA (рис.7.1). Вероятно, на фоне заметного повышения температуры атмосферы только часть этих аномалий отражает изменения притока атлантических вод в Баренцево море. Поэтому в период от 1980 и по 1997 г. индекс N(Y) рассчитывался по величинам TA, как и в более ранних расчетах (Булгакова, 2004), а начиная с 1998 – по величинам TAlin. –отклонениям от линейного тренда (рис. 7.2).

Рисунок 7.1. TA на Кольском разрезе по месяцам, сглаженные значения по трем точкам и линейный тренд за последние 10 лет, причем x=0 соответствует 1 января 1998 г.

Рисунок 7.2. Индекс N(Y), рассчитанный по ежемесячным TA относительно средних температур за период 1951-2000 гг., а начиная с 1998 г. – по отклонениям относительно линейного тренда TA.

Отклонения от линейного тренда рассчитываются по формуле

TAlin =TA– (0,0096*month +0,0106),

где переменная "month" меняется от 0 до 120, причем 0 на оси абсцисс соответствует январю 1998 г.

Рисунок 7.3. Динамика биомассы нерестового запаса SSB(Y), индекса популяционной плодовитости IndPF(Y) и отношения m(Y)=PF/SSB, рассчитанные на основе прогона IUU-Росс XSA (ICES 2008). Здесь Y - год рождения поколения. Величины SSB(Y) отложены по правой оси ординат.

Характер изменчивости популяционной плодовитости и биомассы нерестового запаса похожи, но их отношение не постоянно (зеленая кривая на рис. 7.3), это подтверждает наше утверждение, что в формуле (7-2) параметр m должен меняться с годами, поскольку популяционная плодовитость отражает не только величину, но и структуру нерестового запаса.

В результате параметризации модели пополнения трески получены четыре набора значений параметров (табл. 7.1). Корреляционные отношения для каждого из вариантов показывают, что модель (7-3) описывает от 70% до 77% изменчивости пополнения.

Следует отметить, что для каждого из прогонов XSA, использованных для параметризации модели пополнения, на динамику оценок пополнения Rmod() мало влияет выбранный интервал лет, различия между соответствующими кривыми почти не заметны (рис.7.4). В частности, для прогнозного периода (Y=2001-2007 гг.) различия между ними лежат в пределах 1-3%, а коэффициент корреляции между этими двумя рядами равен 0,999.

Таблица 7.1

Оценки параметров модели пополнения (7-3) и значения корреляционного отношения для разных наборов входных данных, полученных по XSA

Прогон XSA	IUU-Росс	IUU-Норв
Годовые классы, выбран-ные для параметризации	1984 - 2000	1984 - 2003	1984 -2000	1984 -2003
Количество точек / степеней свободы	17/13	20/16	17/13	20/16
Параметры	Оценки параметров
A	4,627	4,709	4,525	4,599
Can	1,535E-06	1,474E-06	1,582E-06	1,497E-06
ny	0,037	0,040	0,032	0,034
Корреляц. отношение r2	0,77	0,72	0,77	0,70

В случае, когда параметризация модели пополнения проводится на основе прогона IUU-Росс XSA, коэффициенты корреляции между рядами R(Y) и Rmod(Y) для обоих интервалов лет равны 0,87, т. е. достаточно велики. Для прогона IUU-Норв они равны 0,88 и 0,82.