CC BY
8Труды ВНИРО 2014 г. Том 151
Водные биологические ресурсы
УДК 639.2.053.7:639.223.5 (265.53)
Оценка запасов и ОДУ минтая восточной части Охотского моря с использованием данных ИС «Рыболовство»
В. К. Бабаян, Т. И. Булгакова, Д. А. Васильев, А. И. Михайлов, И. Н. Антонов,
Г. С. Моисеенко
Всероссийский научно-исследовательский институт рыбного хозяйства и океанографии (ВНИРО, г. Москва) e-mail: vbabayan@vniro.ru
На примере запасов минтая Theragra chalcogramma (Pallas, 1814) восточной части Охотского моря продемонстрированы возможности современного аналитического подхода к обоснованию общего допустимого улова как основной меры регулирования рыболовства. В качестве альтернативного источника исходной информации использована база данных (БД) отраслевой информационной системы (ИС) «Рыболовство», значения уловов на усилие из которой были преобразованы в индексы численности запаса. Oсновное внимание в работе уделялось методическим вопросам, в том числе: стандартизации наблюденных значений уловов на усилие с применением обобщённых линейных моделей (GLM), обоснованию выбора базовой модели динамики промыслового запаса и правил регулирования промысла, риск-анализу и другими диагностическим процедурам, направленным на повышение эффективности управления водными биологическими ресурсами.
Ключевые слова: математическое моделирование, стандартизация уловов на усилие, концепция MSY, правило регулирования промысла (ПРП), минтай Охотского моря.
Введение Одной из основных проблем, препятствующих полноценной оценке состояния запасов промысловых рыб и обоснованию общего допустимого улова (далее ОДУ), является отсутствие или низкое качество доступной информации, используемой в расчётах. В первую очередь это касается данных учетных съёмок, возрастного состава уловов, темпов весового и линейного роста и др. В ряде случаев дефицит необходимой информации не позволяет использовать такие эффективные методы, как различные версии когортного анализа: АДАПТ-метод (ADAPT), статистический анализ уловов по возрастам (SCAA), анализ виртуальной популяции (VPA), расширенный
анализ выживаемости (XSA), а также более сложные методы. В то же время в распоряжении прогнозиста практически всегда имеются промысловые данные, которые в случае интенсивного промысла (каким является промысел охотоморского минтая) могут с успехом заменить результаты научных учетных съёмок и обеспечить информационную основу для применения продукционных моделей.
Уникальным источником промысловых данных в дальневосточном регионе является отраслевая ИС «Рыболовство», в базе данных которой с 1993 г. накапливается информация не только по величине и видовому составу уловов, но и по продолжительности и пространственному распределению промысловых
операций, а также типам рыболовных судов и орудий лова.
В 2011 г. нами была обоснована принципиальная возможность использования данных ИС «Рыболовство» для продукционного анализа запасов западно-камчатского (восточно-охотоморского) минтая и оценки его допустимого промыслового изъятия. Подготовка входной информации для базовой продукционной модели была выполнена с помощью упрощённого метода стандартизации улова на усилия (методом Галланда), информационные требования которого ограничены исходными («сырыми») данными по уловам минтая и соответствующим промысловым усилиям. Полученные таким образом стандартизированные значения уловов на усилие использовались в расчётах в качестве индексов численности. Однако это очень спорное допущение, поскольку на величину стандартизированных таким образом уловов на единицу промыслового усилия (СРЫЕ — далее и) оказывала влияние не только величина запаса, но и целый ряд внешних факторов, которые могли исказить зависимость этого параметра от численности и биомассы запаса.
Для того чтобы приблизить уловы на усилия к индексам численности, в последние годы стал широко применяться метод стандартизации наблюденных значений СРЫЕ, основанный на обобщенных линейных моделях ^ЬМ). Этот метод в противоположность методу Галланда позволяет учесть (и частично нейтрализовать) влияние основных факторов, искажающих истинную зависимость СРЫЕ от величины запаса. Практическая реализация этой задачи для запасов восточно-охотоморского минтая может быть осуществлена с использованием информационных массивов ИС «Рыболовство», которые содержат необходимые для её решения исходные данные.
Основной целью настоящей работы является исследование возможности применения новой методологии, основанной на данных отраслевой информационной системы «Рыболовство».
Краткое описание метода
СТАНДАРТИЗАЦИИ УЛОВА НА УСИЛИЕ
(метода оценки индексов численности)
Процедура стандартизации улова на усилие (Ы) с помощью обобщённых линейных
моделей (метод GLM) представляет собой современную модификацию многофакторного линейного анализа [McCullagh, №Иег, 1989], характерными особенностями которой являются следующие:
применение метода GLM к стандартизации улова на усилие ограничено случаями, при которых распределение иоЬз принадлежит к экспоненциальному семейству распределений случайной величины;
зависимость математического ожидания улова на усилие ^ от линейной комбинации независимых величин (факторов) задаётся некоторой функцией, называемой функцией связи.
Если линейную комбинацию независимых величин (^1, х2, ..., хп), представляющих собой факторы, влияющие на величину и, обозначить как п, то можно записать:
П = р1 X + Р2 *2 + .... + PN*N, (1)
где Р; — неизвестные параметры, характеризующие степень влияния соответствующих независимых переменных х1 на величину улова на усилие; N — число имеющихся записей судовых суточных донесений (далее ССД).
С учётом принятых обозначений обобщенная линейная модель записывается в виде:
фЫ = П или И = ф -1 (п), (2) где ^ — математическое ожидание улова на усилие, а ф — функция связи.
Применительно к задаче стандартизации улова на усилие формула (1) примет вид:
па1,а2.....аК (щ, *2.....%) = У (0 +
+р1*1 +Р2*2 + .... + PN*N + У"1 + (3)
, а2 , , аК
+У2 +... + У К ,
где У (^ — «фактор года», отражающий временную зависимость индекса численности, а уаР — р-й дискретный фактор (р =1,2,...,^, который может принимать несколько значений (категорий), верхний индекс ар — индекс категории.
Поскольку в правую часть уравнения модели помимо временной зависимости фактора численности (биомассы) запаса У (^ входят и другие независимые факторы, определяющие величину СРЫЕ, это позволяет, оценив коэф-
фициенты (в, yapP ) модели (2), выделить факторы (индексы) численности для каждого года промысла из рассматриваемого периода времени. Обычно при правильном выборе независимых переменных степень точности оценки индекса численности растёт с увеличением количества независимых переменных. При этом уменьшается смещение в оценках индекса численности запаса, но возрастает дисперсия этого индекса [Maunder, Punt, 2004].
Анализ и предварительный отбор входных данных
Анализ информации из БД ИС «Рыболовство». Судовое суточное донесение содержит информацию о результатах работы судна за прошедшие отчётные сутки, которые исчисляются по судовому времени, и заполняется на основании документированных данных промыслового, технологического и судового журналов по состоянию на 24:00 по судовому времени. В состав ССД входят следующие показатели, имеющие отношение к промысловой деятельности: код судна в соответствии с реестром судов ИС «Рыболовство»; отчётная дата; координаты; глубина траления; промысловая подзона; объект промысла; суточный улов объекта промысла; орудие лова; суммарный улов за сутки; количество промысловых операций; время, затраченное на промысловые операции; средняя глубина промысловых операций за отчётные сутки.
Из первичных данных, содержащихся в ССД, можно получить значения некоторых вторичных параметров, например, улова на промысловую операцию, улова на час промысловой операции или улов на судосутки промысла. Информация из реестра судов позволяет по коду судна провести разделение промысловых судов по длине корпуса на классы: крупные, средние и малые промысловые суда. Кроме того, возможно деление судов внутри классов на группы по мощности главного двигателя.
В ходе анализа данных ИС «Рыболовство» должны быть решены две задачи: выбор наиболее репрезентативного («стандартного») промыслового комплекса (системы «судно — орудие лова») и обоснование размерности улова на усилие для этого комплекса. Выбранный промысловый комплекс должен отвечать
трём условиям: иметь достаточно равномерное временное (по годам) и пространственное (по районам промысла) распределение, а его доля в общем улове должна быть преобладающей.
Решение второй задачи во многом определяется качеством доступных промысловых данных. На практике при анализе тралового промысла предпочтение часто отдается размерности «улов (в тоннах) на час траления». Однако данные ИС «Рыболовство» по вос-точно-охотоморскому минтаю не позволяют корректно оценить продолжительность тралений для всех зарегистрированных промысловых операций. В этом случае более подходящей размерностью будет «улов (в тоннах) на судосутки промысла». В дальнейшем среднегодовые значения «стандартного» промыслового комплекса с этой размерностью будут преобразованы в индексы численности восточ-но-охотоморского минтая.
Формирование перечня непрерывных и категориальных независимых переменных (факторов-предикторов). По своей природе улов на усилие — это скорее показатель производительности лова, чем показатель численности. В общем случае величина улова на усилие, помимо численности запаса, зависит от двух побочных групп факторов: факторов внешней среды, определяющих поведение особей (плотность пространственного распределения, скорость плавания и др.), и производственных факторов (эффективность орудий и способов лова, профессионализм экипажей промысловых судов, ёмкость трюма и т.д.). Из имеющегося массива первичной информации по ССД отбираются только те характеристики промысла, которые потенциально могут оказывать влияние на величину улова на усилие. При выборе независимых переменных за основу принимается первичная информация из судовых суточных донесений, хранящаяся в базе данных ИС «Рыболовство». Независимые переменные могут представлять собой непрерывные (глубина, координаты траления и др.) или дискретные величины (год, месяц, подзона и др.). При необходимости значения непрерывной величины можно разделить на несколько уровней (категорий), сделав ее дискретной (категориальной). Например, такой
непрерывный фактор, как «глубина траления» легко преобразуется в дискретную переменную путём стратификации значений глубины. Для этого достаточно разделить весь наблюдаемый диапазон глубин на несколько слоёв, присвоить каждому слою определенный индекс и соотнести фактические данные по глубине всех тралений с условными слоями.
Предварительный отбор функций связи для обобщённой линейной модели. Выбор функции связи (ф) модели (2) зависит от закона распределения Ы. В случае обобщённой линейной модели допустимый набор распределений ограничен экспоненциальным семейством (см. Приложение). К экспоненциальному семейству распределений относятся, кроме нормального, некоторые другие распределения, например, обратное нормальное, гамма-, экспоненциальное, а также отдельные дискретные распределения, как например, пуассоновское и биномиальное. Выбор непрерывного или дискретного распределения определяется размерностью зависимой переменной. В нашем случае, где зависимая переменная — улов на единицу промыслового усилия — измерялась в единицах биомассы на промысловое усилие, её распределение с необходимостью было выбрано непрерывным. Но если бы Ы измерялось в численности пойманных особей на единицу промыслового усилия, то при небольших значениях улова на усилие правильнее было бы использовать дискретное распределение.
Выбор потенциально приемлемых функций связи определяется законом распределения наблюденных значений улова на усилие. Например, нормальному распределению и гамма-распределению соответствуют следующие функции связи: тождественная ф(ц) = Ц
—1 Т1
и обратная ф(ц) = Ц . Также существуют распределения, соответствующие логарифмической ф(ц) = 1п ц, степенной (как правило, с показателем 1/2) и логистической функциям. По результатам исследования, представленного в данной работе, наилучшим описанием рассматриваемой базы данных, оказалась обобщённая линейная модель с использованием обратной функции в качестве функции связи.
Стандартизация среднегодовых значений уловов на усилие методом GLM
Выбор «наилучшей» обобщённой линейной модели улова на усилие. Процедура выбора наилучшей модели улова на усилие была построена на серии машинных экспериментов, которые сводились к последовательному тестированию вариантов модели (3), полученных путём различных комбинаций всех предварительно отобранных наборов независимых параметров и версий функций связи. Целью экспериментов являлась оценка адекватности тестируемых вариантов модели исходным значениям уловов на усилие по результатам анализа дисперсий (ANOVA) и сопоставлению значений критерия Акайке, «штрафующим» за избыток неизвестных параметров. Тестирование моделей выполнялось с помощью программных модулей среды Я «ANOVA» и «А1С».
Расчёт стандартизированных значений улова на усилие. Для стандартизации уловов на усилие в отобранную модель следует ввести наиболее представительные значения независимых переменных, которые были найдены путём построения и анализа гистограмм значений каждой из независимых переменных. Чаще всего искомая величина находится как значение данной переменной с наибольшей частотой встречаемости или по наибольшему вкладу в общий улов. В нашем случае оба этих правила выбора приводили к одинаковому результату. В некоторых других случаях в качестве репрезентативного значения предиктора (независимой переменной) было бы оправдано использовать его медиану. Таким образом, оценка стандартизированного индекса биомассы (численности) запаса осуществляется по следующей формуле:
иа () = У(0 + р1 X* +р2*2* + .... (4)
+£n•4 + у*+у* +...+ук
где У(0 — фактор года, х* — «стандартное» (наиболее представительное) значение для всего периода наблюдений непрерывного
фактора, а ур — «стандартное» для всего периода наблюдений значение категориального
фактора, равное величине наиболее представительной категории этого фактора.
Иными словами, стандартизированное значение улова на усилие — это значение линейного предиктора, в котором все переменные, кроме фактора года, постоянны для всех лет рассматриваемого периода времени. В этом и заключается суть стандартизации, выполненной методом GLM.
Эта процедура выполняется для каждого года, по которому имеется необходимая для расчётов информация. Все последующие вычисления осуществляются на основе стандартизированных значений уловов на единицу промыслового усилия как индекса величины запаса. Поэтому для простоты, начиная со следующего раздела, стандартизированные значения уловов на усилие будут обозначаться без индекса <^Ь>, т.е. просто и.
оценка состояния запаса и ВЕЛИЧИНЫ оду (модельный подход)
Оценка параметров продукционной модели. Полученные с помощью метода GLM стандартизированные значения уловов на усилие являются показателями величины запаса в гораздо большей степени, чем исходные «сырые» данные. Это позволяет обоснованно использовать стандартизированные и и соответствующие им величины промыслового усилия Е в качестве входной информации для широкого набора продукционных моделей, на которых основана вся последующая процедура обоснования общего допустимого улова.
В общем виде продукционная модель записывается с помощью двух уравнений, первое из которых выражает динамику запаса, а второе устанавливает связь между биомассой запаса и ее показателем, стандартизированным уловом на усилие I [ВаЬауап, Kizner, 1988]:
В.1+1 = В((1 + гС(В1)) - С;
и. + и.
в =
(5)
2q
где В — биомасса промысловой части запаса в начале года; и — стандартизированный улов на единицу усилия; С — улов за год; С (В) —
функция популяционного роста; г — мгновенный коэффициент весового роста в отсутствии плотностной регуляции (внутренняя скорость роста); q — коэффициент улавливаемости; I — индекс года промысла.
Годовой прирост биомассы запаса описывается функцией популяционного роста С(В). В общем случае С(В) является произвольной монотонно убывающей функцией, имеющей один корень, соответствующий ёмкости среды (величине девственной биомассы запаса) К. При этом наиболее распространено использование двух простых частных случаев:
линейной функции, приводящей к логистическому закону роста С(В) = 1 — В / К ;
логарифмической функции, приводящей к экспоненциальному закону роста
С(В) = -1п(В / К).
Выбор функции роста определяет выбор продукционной модели: в первом случае — это динамический аналог модели Шефера, во втором — модели Фокса.
Выбор типа продукционной модели для использования в дальнейших расчётах можно выполнить двумя способами:
1. По характеру зависимости и(Е). Если фактические данные по усилиям и уловам на усилия лучше аппроксимируются линейной функцией, предпочтение следует отдать аналогу модели Шефера, если экспоненциальной — аналогу модели Фокса.
2. По результатам сравнительного анализа остатков ДЦ (разностей между стандартизированными (и) и теоретическими значениями ( [/ ) улова на усилие, полученными с помощью сравниваемых моделей. Выбор делается в пользу модели, которая дает меньшее суммарное квадратичное отклонение от стандартизированных уловов на усилие.
Оценка неизвестных параметров модели (г, q, К) осуществлялась методом наибольшего правдоподобия путём минимизации предварительно выбранной целевой функции (Ь). Наиболее распространёнными в практике рыбохозяйственных исследований целевыми функциями являются квадратичная, логарифмическая и медианная функции:
N
L=е и - и )2
I=0
N
L = е (1п и. - 1п и. )2
L = med(
и. - и.
где и. — модельное значение улова на усилие в год г; Ц — стандартизированное значение улова на усилие в год г; п — число лет наблюдения.
После нахождения параметров модели выполняется диагностика полученных оценок. С этой целью строятся графики поверхности ошибок в двухпараметрической (плоскостной) проекции и по ним определяется степень устойчивости решений. При построении этих графиков величина одного из параметров модели (г, q, К) фиксируется на оптимальном для этого параметра уровне.
В нашем случае результаты сопоставления всего комплекса расчётных характеристик свидетельствуют в пользу выбора модели Шефера.
Оценка максимального устойчивого улова (ЫЗУ). Сделав в уравнении продукционной модели (5) подстановку и = qB и произведя необходимые преобразования, получим более удобную модификацию модели в виде одного уравнения:
щ -1 )С
и.+. = и. . + г(и. +
г+1 г-1 4 г
Ж и. + и. Л
г-1
2q
2qC..
(6)
В модели (6) состояние запаса характеризуется величиной стандартизированного значения улова на усилие (индексом численности), а интенсивность промыслового воздействия на запас — соответствующим значением промыслового усилия Е. Подставляя в уравнение (6) предварительно преобразованные функции по-пуляционного роста, получим соответствующие динамические версии моделей Шефера и Фокса, выраженные через промысловые характеристики: улов на единицу промыслового усилия и промысловое усилие. Анализ этих моделей в равновесном (статическом) режиме позволяет определить максимальный устойчивый улов
(MSY) и целевые биологические ориентиры по интенсивности промысла (Е^ = EмsY) и биомассе (и1е = UмsY), которые в дальнейшем будут использованы для идентификации правила регулирования промысла (ПРП) на основе концепции MSY. На практике биологические целевые ориентиры обычно выражаются в более привычных величинах: BмsY и FмSY.
Для
перехода к этим величинам используются несложные преобразования: и = qB и Е = . / q.
Оптимизация правила регулирования промысла (ПРП). Для реализации стратегии управления, основанной на регулировании величины годового вылова с учётом концепции MSY, было использовано нелинейное (логистическое) правило регулирования промысла [Бабаян, 2004]. Основное достоинство такого ПРП заключается в том, что, будучи двухзональным, оно не нуждается в граничном ориентире по биомассе, но в то же время обеспечивает повышенную степень защищённости запаса в области низких значений его биомассы. Этим данное правило выгодно отличается от классической 3-зональной схемы предосто-рожного подхода.
Рекомендуемый вариант правила регулирования промысла задается с помощью кусочно-гладкой функции следующим образом: если 0 < В; < 0,5В1г, .гес = (0,5х.^) / (0,5 х
х В1г) а х Вг а; если 0,5 В1г < Вг < В1г, = Гг - (0,5 х
х .1г) / (0,5 х В1г) а х (В1г - В) а; если В; > В1г, .Гес = (7)
где В1 — текущее значение биомассы запаса в год г; Вг и — целевые ориентиры по биомассе и промысловой смертности соответственно (Вг = ВМЗУ + ОВ Ъ; = Гмзу - Ов О, где Ов — стандартная ошибка оценки Bмsy, 1в1 — коэффициент Стьюдента для заданной доверительной вероятности (обычно 90 или 95%); ¥гес г — рекомендуемый уровень промысловой смертности для В = В; а — коэффициент формы.
Если управление направлено на максимизацию среднемноголетнего вылова, то в качестве целевых ориентиров используются значения биомассы и промысловой смертности, соответствующие теоретической величине максимального устойчивого улова. Однако, принимая во внимание неизбежные статистические погреш-
0
Free, год Ftr
Free, i
0 В,
Bi
Br
Рис. 1. Правило регулирования промысла при разных коэффициентах формы (а)
ности в оценках этих параметров, при идентификации ПРП используются более щадящие биологические ориентиры:
Btr = B
■ + Ou L. F„ = F
'МЪУ 1 иВ ^st■, 1 & 1 МвУ °В tst,
где ав — стандартная ошибка оценки Вм$у, ^ — коэффициент Стьюдента для заданной доверительной вероятности (обычно 90 или 95%).
Оптимизация правила регулирования выполняется методом стохастического моделирования динамики запаса (методом Монте-Карло) при различных значениях коэффициента формы "а" определённой выше логистической функции. В качестве критериев оптимизации (эффективности) принятой стратегии управления рекомендуется использовать вероятностные величины, например:
Р(В) — вероятность (риск), что прогнозная биомасса на заданный год будет меньше минимальной биомассы запаса за весь период наблюдений,
Р(С) — вероятность (риск), что средний улов за прогнозный период окажется ниже средней величины улова за такой же по продолжительности период, предшествующий прогнозному.
Для вероятностных критериев эффективности априори устанавливается приемлемый уровень риска, допустимая величина риска по критерию Р(В) обычно не превышает 15%, по критерию Р(С) она составляет не более 25%. Окончательно допустимые уровни риска устанавливаются после консультаций с представителями промышленности.
Оценка рекомендуемой величины ОДУ (модельный подход. Модельная оценка ОДУ+ осуществляется с использованием трехшаговой процедуры:
с помощью динамическои продукционной модели (6) с заданной заблаговременностью k прогнозируется состояние запаса; состояние запаса может быть охарактеризовано в терминах биомассы (B +k) или улова на усилие
(U+k);
по найденному состоянию запаса в год (i+k) с помощью правила регулирования промысла (7) оценивается рекомендуемая величина интенсивности промысла (соответственно
Frec, i+k или Erec, i+k);
рассчитывается величина ОДУ на (i + ^-й год:
ОДУ + = Frec, i+k Bi+k
или
ОДУ^+к = Erec, i+k Ui+k
Изложенная в разделе 4 комплексная процедура оценки запаса и обоснования ОДУ выполнена с помощью специализированного пакета прикладных программ Combi [Бабаян и др., 2011].
Оценка рекомендуемой величины ОДУ
(НЕМОДЕЛЬНЫЙ ПОДХОД)
Если период наблюдений недостаточно продолжителен, а динамика запаса сравнительно стабильна и не отличается значительными колебаниями, то оценить параметры продукционной модели с необходимой точностью и воспользоваться описанной в разделе 4 процедурой не всегда удаётся. В этом случае для оценки ОДУ можно использовать немодельные (эмпирические) методы, которые основаны на учёте тенденции в динамике величины запаса. Тенденция оценивается с помощью количественного анализа изменений в индексах численности за период, непосредственно предшествующий прогнозному. Количественно тенденция оценивается по величине (и знаку) тангенса наклона графика линейной регрессии, построенного в координатах «индекс численности» — «год промысла». Ниже приведены два немодельных метода, которые могут быть рекомендованы для практического применения.
Наиболее простой из этих методов [Bentley, 2002] основан на коррекции предыдущего значения общего допустимого улова (ОД^ _i) путём учёта градиента изменения величины запаса (индекса численности) относительно
тыс. т
некоторого интервала времени, предшествующего прогнозному году г + 1:
ОДУ ,+1 = (I(п) / I, )ОДУ,,
где 1^п) — среднее значение индекса численности за период с (г — п — 1)-го по г-й год включительно. Период усреднения не должен быть слишком большим (так, авторы метода приняли п = 3).
1 г'
I(п) =1 е I 1 п ^ У
У=1-п+1
Второй рекомендуемый метод более полно и наглядно учитывает наличие и характеристики тренда в динамике запаса [ВиИетойЬ,
Geromont, 1997].
Базовое уравнение метода:
ОДУ;+1 = ОДУг (1 + Хв),
где г — индекс года промысла; в — мера крутизны тренда индекса величины запаса (тангенс наклона графика линейной регрессии логарифмов индексов численности за последние несколько лет); X — корректирующий коэффициент, устанавливающий уровень чувствительности регулирования.
Согласно практике применения этого метода в Организации по рыболовству в Северо-Западной Атлантике (НАФО) [NAFO, 2010] при управлении запасами гренландского палтуса в зависимости от направления тренда коэффициент X может принимать значения 1,00 (если запас растёт, в > 0) и 1,25 или 2,00 (если запас стабилен или снижается, в < 0). Временной интервал, на котором оценивается коэффициент в, принят равным 5 годам.
При применении немодельных методов отдельного решения требует задача выбора стартового значения улова. В случае, если в предшествующие годы доступная информация позволяла осуществлять полное аналитическое обоснование ОДУ, выбор должен быть однозначным — последняя оценка ОДУ. Если на протяжении последних лет оценка ОДУ осуществлялась эмпирическими методами, но при этом заметного снижения фактического вылова не наблюдалось, то в качестве стартового значения улова следует принять среднее значение ОДУ за несколько предшествующих лет. Если в предшествующий период анали-
тическая оценка ОДУ не осуществлялась, то в качестве стартовой величины улова можно использовать величину фактического улова в предпрогнозный год (если динамика вылова достаточно стабильна) либо среднюю величину фактического вылова за предшествующий период (если наблюдаются значительные среднегодовые колебания уловов).
результаты расчётов и их ОБСУЖДЕНИЕ
Исходные данные. В качестве исходных данных использовалась информация судовых суточных донесений из базы данных ИС «Рыболовство» за март 1998—2010 гг., т.е. за 14-летний период наблюдений. На март приходится пик добычи восточно-охотоморско-го минтая, в связи с чем промысловые данные в этом месяце наиболее полно характеризуют плотность и численность облавливаемого запаса. Исходная информация была ограничена данными по разноглубинным тралам, поскольку этим типом орудия лова на протяжении рассматриваемого периода вылавливалось до 83% от общего улова, что позволяет считать разноглубинные тралы репрезентативным орудием лова для промысла восточно-охотоморского минтая. При этом учитывалась промысловая информация только по тем тралениям, в которых минтай составлял не менее 50% улова. Фрагмент исходных данных приведен в табл. 1.
Формирование перечня непрерывных и дискретных переменных (факторов-предикторов). В качестве зависимой переменной рассматривается улов на единицу промыслового усилия с размерностью «улов (т) на судо-сутки промысла».
В качестве независимых переменных (предикторов) выбраны следующие факторы: класс судна, год, район промысла, орудие лова, глубина траления, мощность главного двигателя промыслового судна.
Согласно официальному районированию восточной части Охотского моря фактор «Район промысла» (Area) включает 2 категории: 272 — Камчатско-Курильская подзона, 274 — Западно-Камчатская подзона.
Факторы «Класс судна» и «Мощность судна» были объединены в один фактор «Class/ Power», учитывающий обе эти характеристики.
Таблица 1. Характеристики промысловых операций (тралений) российских судов (фрагмент БД ИС «Рыболовство»). В первой строке указаны названия факторов в GLM
Year Area CPUE Depth Type_code Class_code FG Power
1998 272_25_130_171_8 tral251-500_588
1998 272_30_420_128_7_tral751_5700
1998 272_123_108_171_8 tral251-500_589
1998 272_5,14_110_170_8 tral501-750_971
1998 274 127,37 500_168_8 tral501-750 1380
1998 272 123,6_55_168_8 tral501-750 1380
1998 274_141,5 500_168_8 tral251-500 1380
1998 274_92,6_430_168_8 tral501-750 1380
1998 274 105,5 530 168 8 tral501-750 2425
Фактор «Орудие лова» был разделен с учётом длины гужа (I) разноглубинного трала на четыре категории: I < 250 м, 250 < I < 500 м, 500 < I < 750 м и I > 750 м.
Фактор «Глубина траления» стратифицирован на 11 категорий (табл. 2)
Таблица 2. Категории фактора «Глубина траления» (D)
1. D < 100 м 7. 600 < D < 700 м
2. 100 < D < 200 м 8. 700 < D < 800 м
3. 200 < D < 300 м 9. 800 < D < 900 м
4. 300 < D < 400 м 10 . 900 < D < 1000 м
5. 40 < D < 500 м 11 . D > 1000 м
6. 500 < D < 600 м
Отбор функций связи для обобщённой линейной модели улова на усилие. Поскольку строгого правила определения закона распределения случайной величины не существует, по аналогии с выбором закона распределения улова на усилие, приведенного в ряде опубликованных работ, принимаем гипотезу о том, что в нашем случае распределение и относится к семейству экспоненциальных распределений. Наиболее вероятными распределениями в нашем случае будут гамма-распределение и обратное гауссовское распределение, поскольку они не симметричны. Для указанных распределений функция связи может быть выбрана из следующих функций: тождественной, логарифмической и обратной. Окончательный вы-
бор функции связи осуществляется по результатам анализа стохастического моделирования версий обобщённой линейной модели улова на усилие с различными наборами функций связи и независимых факторов.
стандартизация среднегодовых значений уловов на усилие методом GLM. Выбор «наилучшей» модели CPUE. В машинных экспериментах использовался следующий вид уравнения для математического ожидания наблюдавшегося ранее улова на усилие вида:
EUobs = j-1 (Year + Area + + FG + Depth + Class/Power),
где Year —коэффициент фактора «Год»; Area — коэффициент фактора «Район промысла»; FG — коэффициент фактора «Орудие лова», Depth — коэффициент фактора «Глубина траления»; Class/Power — коэффициент фактора «Класс и мощность судна».
Для выбора «наилучшей» модели улова на усилие были протестированы версии модели с предварительно выбранными функциями связи (тождественная, логарифмическая, обратная) и различным набором независимых переменных (факторов). Отметим, что в случае категориальных факторов в модели были представлены все категории в качестве независимых параметров. По результатам машинных экспериментов были проведены статистические тесты: анализ дисперсии (ANOVA) и расчёт информационного критерия Акайке (AIC). Тесты показали, что «наилучшей» ока-
залась модель с обратной функцией («inverse» в пакете R) в качестве функции связи и полным набором учитываемых факторов.
Расчёт индексов численности. Для расчёта годовых индексов численности (стандартизированных уловов на усилие) использовалось следующее уравнение:
U J = ф-1 (Year + Areast+ FGst + + Depthst + ClassPowerst).
При выборе стандартных категорий дискретных факторов для оценки стандартизированных уловов на усилие учитывались следующие показатели:
величина вылова, соответствующая каждой категории рассматриваемого дискретного фактора;
частота встречаемости каждой категории рассматриваемого фактора в статистических районах восточной части Охотского моря 1.
На основании приведённых выше критериев в качестве стандартной категории фактора «Орудия лова» были выбраны разноглубинные тралы, у которых длина гужа находится в диапазоне 500—750 метров. Количество этих тралов превышает 46% от общего числа тралов, использованных на промысле вос-точно-охотоморского минтая за весь период наблюдений.
В качестве стандартной категории фактора «Класс и мощность судна» выбрано среднетоннажное судно с главным двигателем мощностью 900 — 1000 л. с. Данная категория промысловых судов преобладала во
все годы промысла и составляла 15—38% от общего числа судов. Эта категория внесла наибольший вклад в общий вылов за весь рассматриваемый период. Кроме того, промысловые суда этой категории численно преобладали в обоих промысловых районах. По тем же критериям были выбраны стандартные категории факторов «Глубина траления» (200—300 м) и «Район промысла» (274). Результаты стандартизации уловов на усилие приведены в табл. 3 и на рис. 2.
45
40
ш ? 35
а
и
30
и
25-
I П
III
1998 2000 2002 2004 2006 2008 2010
Год
Рис. 2. Результаты GLM-оценки индекса численности и доверительных интервалов (при p = 0,95)
Оценка состояния запаса и величины ОДУ (модельный подход). Оценка параметров продукционной модели. Все последующие расчёты были выполнены при помощи пакета прикладных программ Combi.
Таблица 3. Оценка индекса численности и промыслового усилия на 1998—2010 гг.
Год Индекс численности, т/судосутки промысла Промысловое усилие, судосутки промысла, 103 Год Индекс численности, т/судосутки промысла Промысловое усилие, судосутки промысла, 103
1998 40,56 20,81 2005 30,89 8,59
1999 34,81 14,62 2006 32,25 9,81
2000 38,21 8,98 2007 34,58 8,52
2001 33,07 10,14 2008 34,26 13,25
2002 30,35 7,95 2009 33,93 15,37
2003 35,2 8,00 2010 37,71 19,19
2004 30,9 5,87
1 Для выбора стандартной категории Агеа^ учитывался только первый показатель.
Исходная информация для использования продукционной модели такова: расчётные годовые значения стандартизированного улова на усилие и за период 1998—2010 гг. и соответствующие значения промысловых усилий Е51 = С / и^, где С^ — общий годовой вылов минтая в восточной части Охотского моря всеми типами промысловых судов и орудий лова.
На основании статистических критериев в качестве базовой модели выбрана динамическая версия продукционной модели Шефера и квадратичная целевая функция. При оценке параметров продукционной модели максимально допустимое значение одного из параметров — ёмкости среды (К), было ограничено величиной 3400 тыс. т [Шунтов, Дулепова, 1997]. С учётом данного ограничения были получены следующие оценки:
внутренний (популяционный) темп роста г = 0,61;
коэффициент улавливаемости q = 0,0124; ёмкость среды К = 3391 тыс. т. Оценка максимального устойчивого улова (ЫЗУ). Максимальный устойчивый улов рассчитан с использованием статической версии динамической продукционной модели и составил 499,7 тыс. т. Соответствующие значения промысловой смертности и биомассы запаса (В) равны: FмsY = 0,228 год-1, = 1699,6 тыс. т.
Идентификация правила регулирования промысла (ПРП). Согласно расчётам биомасса запаса в терминальном 2010 г. равна 2891 тыс. т, что почти в два раза превышает ВМ5у. В этом случае в соответствии с зональным правилом регулирования промысла (п. 4.3, (7)) и концепцией МБУ рекомендуемая промысловая смертность будет совпадать с целевым ориентиром: Ргес = = = Емзу. Для этого случая целевое значение промысловой смертности было рассчитано как нижняя граница доверительного интервала оценки РШу при доверительной вероятности Р = 90%:
Р^мзу - = 0,199 год-1, где aF — стандартная ошибка оценки Рм$у, ,, — коэффициент Стьюдента при доверительной вероятности Р = 90%.
Оценка рекомендуемой величины ОДУ. Общий допустимый улов рассчитывался по формуле: ОДУ; = Fмsу В. При прогнозировании В!, начиная с г = 2012, улов за предшествующий год принимался равным оценке ОДУ на этот год. Результаты оценки ОДУ и математического ожидания биомассы запаса представлены в табл. 4 и на рис. 3.
На рис. 3 представлены ретроспективная и прогнозируемая динамика биомассы и ОДУ в условиях предложенного управления. Из графика видно, что предложенная схема регулирования промысла приводит к стабилизации биомассы на уровне, близком к оптимально-
Таблица 4. Оценка ОДУ и биомассы запаса на 2011—2020 гг.
„ Биомасса запаса, „ , „ Биомасса запаса, „ ,
1од ОДУ, тыс. т 1од ОДУ, тыс. т
тыс. т тыс. т
2011 2525 742 2021 1729 508
2012 2247 661 2022 1722 506
2013 2077 611 2023 1716 505
2014 1966 578 2024 1712 503
2015 1890 556 2025 1709 502
2016 1838 540 2026 1706 502
2017 1800 529 2027 1705 501
2018 1774 521 2028 1703 501
2019 1754 516 2029 1702 500
2020 1740 512 2030 1702 500
1995 2000 2005 2010 2015 2020 2025
Год
Рис. 3. Динамика промысловой биомассы запаса минтая FSB (с 2011 г.— прогноз), полученная по продукционной модели, общей биомассы по когортной модели Syntesis TSB [Антонов, 2011], официальные уловы и ОДУ (начиная с 2010 г.)
му в смысле обеспечения наибольшего улова в долгосрочной перспективе. На прогнозном участке откладываются не сами реальные величины биомассы и ОДУ, а только их математические ожидания, поэтому флуктуации отсутствуют, поскольку их принципиально невозможно предсказать заранее, можно только оценить их дисперсию.
Из рисунка 3 видно, что управление согласно предложенному правилу регулирования позволило бы в 2010—2012 гг. получить более высокий улов по сравнению с фактическим. Оценки биомассы запаса, полученные по продукционной модели, несколько ниже оценок, полученных по модели Synthesis. Различия объясняются тем, что в первом случае это оценка промыслового запаса, во втором — общего запаса.
Оценка ОДУ (немодельный подход).
Для оценки ОДУ минтая на 2011—2012 гг. использовался также немодельный метод (раздел 5). Для оценки в было построено уравнение линейного тренда логарифмов индексов численности за последние пять лет наблюдения (2006—2010 гг.). Наклон тренда (в) составил 0,029. Оценка ОДУ2011 = 744,9 тыс. т. получена по формуле:
ОДУ2011 = С2010 (1 + 0,029),
где С2010 — фактический вылов за 2010 г., равный 723,6 тыс. т.
Индекс численности на 2011 г. рассчитан путём построения линейной регрессии U на ОДУ. Коэффициент в был пересчитан на период 2007—2011 гг. и составил 0,024. Оценка ОДУ2012 = 762,9 тыс. т. получена по формуле:
ОДУ2012 = ОДУ2011 (1 + 0,024),
где ОДУ2011 — оценка ОДУ на 2011 г., равная 744,9 тыс. т.
Выводы
Предположенная процедура обработки промысловой информации отраслевой базы данных на основе современной методологии обобщённых линейных моделей позволила получить ряд репрезентативных индексов численности и на его основе провести полный анализ системы «запас — промысел», обеспечивающий в конечном счёте оценку ориентиров управления и идентификацию правила регулирования промысла. Расчёты показали, что текущее состояние запаса восточно-охо-томорского минтая и его прогнозируемое на ближайшую перспективу состояние в условиях предложенного управления можно охарактеризовать как устойчивое, которое позволяет поддерживать промысел на уровне, близком к теоретическому максимуму годовой продуктивности запаса (MSY).
Работа поддержана грантом РНФ 1411-00687.
ЛИТЕРАТУРА
Антонов Н. П. 2011. Промысловые рыбы Камчатского края: биология, запасы, промысел. М.: Изд-во
ВНИРО. 244 с.
Бабаян В.К. 2004. Альтернативные методы оценки рекомендуемой интенсивности промысла при расчёте ОДУ // Рыбное хозяйство. Т. 4. С. 23—25. Бабаян В. К., Антонов И. Н., Михайлов А. И. 2011. Программный комплекс «СотЫ». Свидетельство об официальной регистрации программ для ЭВМ № 2011615622. Реестр программ для ЭВМ. Бабаян В. К., Варкентин А.И., Васильев Д. А., Сергеева Н. П. 2006. Методические особенности обоснования ОДУ минтая в условиях неопределённости // Методические аспекты исследований рыб морей Дальнего Востока. Труды ВНИРО. Т. 146. М: Изд-во ВНИРО. С. 13-36. Шунтов В. П., Дулепова Е. П. 1997. Современный статус, био- и рыбопродуктивность экосистемы Охотского моря // Комплексные исследования экосистемы Охотского моря. М.: Изд-во ВНИРО.
С. 248-261.
Babayan V, Kizner Z. 1988. Dynamic Models for TAC Assessment: Logic, Potentialities, Development //
Colln. Scient. Pap. Int. Commn. SE. Atl. Fish. 15 (I). Р. 69-83.
Babayan V., Antonov I., Bulgakova T., Rolskiy A., Tretyakov I., Vasilyev D. 2012. The Assessment of Stock Status and Commercial Potential of Redfish in the Irminger Sea // NEAFC, АМ 2012-37. 20 р. McCullagh P., Nelder J. 1989. Generalized Linear Models. Second Edition. Boca Raton: Chapman and Hall /
CRC. IBSN 0-412-31760-5. Maunder M. N, Punt A. E. 2004. Standardizing Catch and Effort Data: A Review of Recent Approaches // Fisheries Research. V. 70. Р. 141-159. Butterworth D., Geromont H. 1997. Evaluation of a Range of Possible Simple Interim Management Procedures of the Namibian Hake Fishery // Report to the Ministry of Fisheries and Marine Resources. Namibia. 28 p. Kell L. T., Mosqueira I., Grosjean P., Fromentin J-M., Garcia D., Hillary R., Jardim E., Mardle S., Pastoors M. A, Poos J. J, Scott F., Scott R. D. 2007. FLR: An Open-Source Framework for the Evaluation and Development of Management Strategies // J. of marine sci. Vol. 64. № 4. P. 640-646. NAFO. 2010. Report of the Working Group on Greenland Halibut Management Strategy Evaluation (WGMSE) / Serial No. 5866 NAFO/FC Doc. 10/30. 64 p. Bentley N., Breen P. A, Starr P. J., Sykes D. R. 2002. Development and Evaluation of Decision Rules for Management of New Zealand Rock Lobster Fisheries // New Zealand Fisheries Assessment Report.
Приложение свойства экспоненциального семейства. Стандартизация уловов на усилие с использованием метода GLM (если рассматривать её место в рамках процедуры MSE [Kell et al., 2007]) есть не что иное, как настройка модели наблюдения. Модель наблюдения — это статистическая гипотеза о распределении зависимой переменной. Обобщённая линейная модель предписывает рассматривать распределения наблюдаемой величины только из экспоненциального семейства, общий вид которого задаётся следующим образом [McCullagh, Nelder,1989]:
exp(u|- L(u) - G(r|)). (A.1)
Величина n называется предиктором модели. Все статистические моменты наблюдаемой величины будут функциями предиктора. Статистическая модель называется обобщённой аддитивной, если предиктор как функция
многих независимых переменных аддитивно факторизуема, т.е. представляет собой сумму функций одной переменной:
= у(,) + № + /2(х2) + ... + !п(хп) + (А.2) +У. + Уг()2) +... + УтЮ.
Независимые переменные могут пробегать как непрерывное, так и дискретное множество значений. В последнем случае такие переменные называются категориальными. В GLM принято предположение о линейной зависимости предиктора от непрерывных факторов:
, х2хп; .г.) =
= у(1) + Ь +Р2 Х2... + РпХп + (А.3) +У1О1 ) + У2 ) + ... + Ут (т).
Коэффициенты при факторах у(,), {Р^},{у (. )}являются неизвестными параметрами, которые необходимо определить при настройке модели. Благодаря факторизации количество наблюдений значительно превышает число оцениваемых параметров, и именно поэтому аппарат статистики может быть применён.
В уравнении (А.1) Ь(и) есть взятая с обратным знаком функция правдоподобия (логарифм плотности вероятности) некоторого, достаточно произвольного распределения, а
ехр(С(ц)) = Техр(иц - Ь(и)^и
есть производящая функция моментов этого распределения. Таким образом, для всякого распределения, обладающего конечной производящей функцией моментов, можно построить содержащее его экспоненциальное семейство.
Экспоненциальное семейство указанного вида обладает следующими замечательными свойствами:
1. Зависимость между математическим ожиданием и линейным предиктором однозначно определяется обратной функция связи:
, . дС(ц) -и , дц
2. Отклонения наблюдений от математического ожидания имеют вид:
3. Дисперсия (в общем случае — матрица ковариаций) является функцией линейного предиктора, а значит, и математического ожидания:
D(u) = Э2ЛЛОД =
д2
М(п) = д пЕ ехр(ги)
г=0
и ж ехр (-
(и - 2)2
а
1 N
2=- е и..
N ^ ¡
Предположим, что:
2=2о+еь
клк
(А.4)
где Xk — факторы, а в к — параметры. Дифференцируя по параметрам функцию правдоподобия:
N
ь ~ е (и -2о -еРЛ )2>
¡=1
k
4. Производящая функция моментов:
Е ехр(Ш) = ехр(С(2 + t) - 6(п)).
Для вычисления любого момента функции распределения достаточно соответствующее число раз продифференцировать в нуле производящую функцию:
5. Уравнение регрессии определяется только функцией связи и имеет следующий вид:
е^) д|=
Это уравнение, в отличие от уравнения линейной регрессии, нелинейно относительно параметров. Уравнение линейной регрессии является частным случаем данного уравнения, когда функция связи линейна.
Перечисленные выше свойства экспоненциального семейства распределений обусловлены тем, что функция связи однозначно задает вид распределения и все его моменты.
Поскольку многофакторная регрессия является частным случаем обобщенных линейных моделей (при линейной функции связи), то наиболее простым примером будет гипотеза о нормальном распределении наблюдаемой величины. В этом случае задача стандартизации сведется к тривиальному усреднению. Действительно,
легко получить хорошо известную систему регрессионных уравнений на параметры модели:
N
дРь = е(и - 2°х -ЕРлх) = °;
I=1 k
N
д ь=е (u¡ -2° ^Хк)=°.
¡=1 k
Другой способ обобщения заключается в отказе от квадратичной функции правдоподобия. Для этого перепишем функцию правдоподобия нормального распределения в виде:
ь - (и - V)2 = 1 и2 - и2- 122.
Теперь сделаем функцию распределения произвольной по отношению к независимой переменной, но определённой по отношению к предиктору:
ь(и, 2) = ь(и) - и2 + G(2).
Функция G(n) определяется из условия нормировки:
| ехр(-Ь(и, 2)) ¿и = 1"2 или, что то же самое:
G(2) = 1п Т ехр(-Ь(и) + и2)^и.
Дифференцирование функции правдоподобия приводит к уравнению:
N
дць = ° Ы е и = д^(2).
¡=1
С другой стороны, с учётом:
Т иехр(-Ь(и) + и2^и
т(2) = дЛОД = -= Еи
| ехр(-Ь(и) + и2)^и
вид функции связи однозначно определён видом распределения.
Для того чтобы завершить построение обобщенной линейной модели, достаточно
к
предположить линейную связь зависимой переменной и параметров в виде (А.4), а также рассматривать переменные и, п и хк не как скалярные, а как векторные величины одинаковой размерности.
Следует отметить, что распределение зависимой переменной в обобщенной линейной модели не является совершенно произвольным. Так, например, функция правдоподобия вида L ~ (и — ц) не принадлежит классу обоб-
щённых линейных моделей, так как содержит члены вида ип2 и и2п. Тем не менее, интерес к использованию обобщённых линейных моделей обусловлен более широким классом распределений моделируемой переменной, поэтому GLM можно успешно применять в том случае, если оценка по многофакторной регрессии несостоятельна, т.е. гипотеза о нормальном распределении зависимой переменной признается статистически недостоверной.
Walleye Pollock of the East Part of the Okhotsk Sea Stock Size and TAC Assessment with the Use of «Rybolovstvo» Data Base
V. Babayan, T. Bulgakova, D. Vasilyev, A. Mikhaylov, I. Antonov, G. Moiseenko
Federal Research Institute of Fisheries and Oceanography (VNIRO, Moscow) e-mail: vbabayan@vniro.ru
The potenciáis of the modern analytical approach to assessment of a total allowable catch (TAC) as the key fisheries management measure were demonstrated on the example of the stock of walleye pollock Theragra chalcogramma (Pallas 1814) in the East of the Okhotsk Sea. The database (DB) of the information system (IS) «Fishery» was used as an alternative source of the initial information for analysis. The observed catches per unit effort ware converted into stock abundance indices. The main emphasis was placed on the following methodological issues: standardization of the observed values of catch per unit effort using generalized linear models (GLM), justifying the choice of the base model of the stock dynamics and the harvest control rule (HCR), risk analysis and other diagnostic procedures aimed to improve the efficiency of the aquatic biological resources management.
Key words: mathematical modeling, CPUE standardization, IS «Rybolovstvo», MSY concept, harvest control rule, Okhotsk sea, walleye pollock.
