17I. DLM: Введение в методы анализа при очень ограниченных данных
Это занятие открывает серию практик о так называемых “немодельных”, как написано в некоторых методиках, подходах, которые в англоязычной литературе принято называть DLM. Это занятие вводное - возьмем cамый простой метод и прогоним весь цикл. Data-Limited Methods (DLM) — это не просто статистические подходы, а скорее философия работы в условиях неопределенности. Если классические методы оценки запасов требуют многолетних рядов данных, детальной информации о возрасте, размерах, промысловом усилии и прочих параметрах, то DLM предлагают прагматические решения для ситуаций, когда единственное, что у вас есть — это данные об уловах. И да, такое бывает чаще, чем кажется.
Представьте себе типичную ситуацию: новый объект промысла, развивающийся рыболовный регион, или просто вид, который никогда не был приоритетом для исследований. Данные об уловах есть — их собирают всегда. А вот биологические параметры, индексы численности, данные научных съемок — это уже из области фантастики. Именно в таких условиях DLM становятся не просто удобным инструментом, а единственным способом хоть как-то оценить состояние запаса.
Среди множества DLM-методов особое место занимает Catch-MSY — подход, который позволяет оценить максимальный устойчивый вылов (MSY) на основе одних только данных об уловах. Метод основан на простой, но гениальной идее: если мы знаем историю промысла, мы можем методом подбора найти такие комбинации параметров роста популяции (r и K), которые объясняют наблюдаемую динамику уловов и при этом биологически реалистичны.
Конечно, у этого метода есть ограничения. Он предполагает, что популяция следует логистической модели роста, что промысел был основным фактором смертности, и что мы можем разумно ограничить диапазоны параметров *r** и K. Но как говорится, в условиях слепоты и одноглазый — король.
Что особенно важно — Catch-MSY не просто дает точечную оценку MSY, но и обеспечивает распределение вероятностей различных значений, что позволяет оценить неопределенность. Это важно для управления рисками.
Переход от оценки к управлению — это отдельная история. Полученные оценки MSY становятся основой для разработки правил управления (Harvest Control Rules), которые определяют, как следует изменять уровень изъятия в зависимости от состояния запаса. А уже эти правила тестируются в рамках Management Strategy Evaluation (MSE) — своеобразного полигона, где мы проверяем, как различные стратегии управления будут работать в условиях неопределенности.
Курьезность ситуации заключается в том, что несмотря на то, что мы используем чрезвычайно простые методы для оценки, процесс управления становится все более сложным. Но это именно тот случай, когда сложность оправдана — ведь на кону устойчивость промысла и сохранение ресурсов.
Ирония судьбы: методы, разработанные для условий недостатка данных, часто оказываются более робастны (устойчивы), чем их сложные собратья, требующие тонны информации. Просто потому что они “осознают” неопределенность, а не делают вид, что ее не существует.
В следующих лекциях мы перейдем от теории к практике и посмотрим, как эти методы работают на реальных данных. Приготовьтесь к приятным сюрпризам.
# =============================================================================# ПРАКТИЧЕСКИЙ АНАЛИЗ ВОДНЫХ БИОРЕСУРСОВ МЕТОДОМ CATCH-MSY# =============================================================================# ------------------------- 1. УСТАНОВКА И ЗАГРУЗКА ПАКЕТОВ ---------------------------# Устанавливаем необходимые пакеты (если еще не установлены)if (!require("ggplot2")) install.packages("ggplot2")
Загрузка требуемого пакета: ggplot2
if (!require("dplyr")) install.packages("dplyr")
Загрузка требуемого пакета: dplyr
Присоединяю пакет: 'dplyr'
Следующие объекты скрыты от 'package:stats':
filter, lag
Следующие объекты скрыты от 'package:base':
intersect, setdiff, setequal, union
if (!require("tidyr")) install.packages("tidyr")
Загрузка требуемого пакета: tidyr
# Загружаем пакетыlibrary(ggplot2)library(dplyr)library(tidyr)# ------------------------- 2. ЗАГРУЗКА ДАННЫХ ---------------------------# Вектор лет наблюденийYear <-2005:2024# Данные по вылову (тыс. тонн)Catch <-c(5, 7, 6, 10, 14, 25, 28, 30, 32, 35, 25, 20, 15, 12, 10, 12, 10, 13, 11, 12)# Создаем график динамики выловаcatch_df <-data.frame(Year = Year, Catch = Catch)ggplot(catch_df, aes(x = Year, y = Catch)) +geom_bar(stat ="identity", fill ="steelblue", alpha =0.7) +geom_line(color ="red", linewidth =1) +labs(title ="Динамика вылова водных биоресурсов (2005-2024)",y ="Вылов, тыс. тонн", x ="Год") +theme_minimal() +theme(plot.title =element_text(hjust =0.5))
# ------------------------- 3. РЕАЛИЗАЦИЯ МЕТОДА CATCH-MSY ---------------------------# Функция для реализации метода Catch-MSYcatch_msy_analysis <-function(catch, n_iter =100000, r_prior =c(0.01, 1.5)) { n_years <-length(catch)# Приоры для параметров r <-runif(n_iter, r_prior[1], r_prior[2]) k <-runif(n_iter, max(catch) *3, max(catch) *50) msy <- r * k /4# MSY = r*K/4 для модели Шефера# Массивы для хранения результатов viable <-rep(FALSE, n_iter) b_msy <- k /2# Bmsy = K/2 для модели Шефера# Проверяем каждую комбинацию параметровfor (i in1:n_iter) { biomass <-numeric(n_years +1) biomass[1] <- k[i] # Начальная биомасса = K# Моделируем динамику биомассыfor (t in1:n_years) { surplus_production <- r[i] * biomass[t] * (1- biomass[t] / k[i]) biomass[t +1] <-max(0.001* k[i], biomass[t] + surplus_production - catch[t]) }# Проверяем условия жизнеспособностиif (all(biomass >0.001* k[i]) &&all(biomass <1.1* k[i]) && biomass[n_years +1] >0.2* k[i] && biomass[n_years +1] <0.8* k[i]) { viable[i] <-TRUE } }# Возвращаем жизнеспособные комбинации параметровlist(r = r[viable],k = k[viable],msy = msy[viable],b_msy = b_msy[viable],depletion = biomass[n_years +1] / k[viable] )}# Применяем метод Catch-MSYset.seed(123) # Для воспроизводимости результатовmsy_results <-catch_msy_analysis(Catch)# ------------------------- 4. АНАЛИЗ РЕЗУЛЬТАТОВ CATCH-MSY ---------------------------# Основные статистикиmsy_summary <-data.frame(Parameter =c("r", "K", "MSY", "Bmsy"),Median =c(median(msy_results$r),median(msy_results$k),median(msy_results$msy),median(msy_results$b_msy) ),Mean =c(mean(msy_results$r),mean(msy_results$k),mean(msy_results$msy),mean(msy_results$b_msy) ),SD =c(sd(msy_results$r),sd(msy_results$k),sd(msy_results$msy),sd(msy_results$b_msy) ))print("Результаты анализа Catch-MSY:")
[1] "Результаты анализа Catch-MSY:"
print(msy_summary)
Parameter Median Mean SD
1 r 0.09671814 0.1409528 0.1330245
2 K 487.67637851 555.3402905 292.9906278
3 MSY 13.08343781 12.7539980 5.7560545
4 Bmsy 243.83818925 277.6701453 146.4953139
# Визуализация распределения параметровpar(mfrow =c(2, 2))hist(msy_results$r, main ="Распределение r", xlab ="r", col ="lightblue")hist(msy_results$k, main ="Распределение K", xlab ="K", col ="lightgreen")hist(msy_results$msy, main ="Распределение MSY", xlab ="MSY", col ="lightcoral")hist(msy_results$depletion, main ="Распределение деплетированности", xlab ="Bcurrent/K", col ="lightyellow")
# Текущий статус запасаcurrent_status <-data.frame(Metric =c("Текущий вылов", "MSY", "Отношение вылова к MSY", "Деплетированность (B/K)"),Value =c(mean(Catch[length(Catch)-2:0]), # Средний вылов за последние 3 годаmedian(msy_results$msy),mean(Catch[length(Catch)-2:0]) /median(msy_results$msy),median(msy_results$depletion) ))print("Текущий статус запаса:")
[1] "Текущий статус запаса:"
print(current_status)
Metric Value
1 Текущий вылов 12.0000000
2 MSY 13.0834378
3 Отношение вылова к MSY 0.9171901
4 Деплетированность (B/K) 3.1261895
# ------------------------- 5. ОПРЕДЕЛЕНИЕ REFERENCE POINTS ---------------------------# Reference points на основе Catch-MSYreference_points <-data.frame(Point =c("MSY", "Bmsy", "Fmsy"),Value =c(median(msy_results$msy),median(msy_results$b_msy),median(msy_results$r) /2# Fmsy = r/2 для модели Шефера ))print("Reference Points:")
[1] "Reference Points:"
print(reference_points)
Point Value
1 MSY 13.08343781
2 Bmsy 243.83818925
3 Fmsy 0.04835907
# ------------------------- 6. ПОСТРОЕНИЕ ПРАВИЛА УПРАВЛЕНИЯ (HCR) ---------------------------# Функция для определения HCR на основе B/Bmsyharvest_control_rule <-function(b_ratio, blim =0.4, btarget =0.8, fmin =0.05, fmax =1.0) {if (b_ratio <= blim) {return(fmin) # Минимальный уровень изъятия при низкой биомассе } elseif (b_ratio >= btarget) {return(fmax) # Максимальный уровень изъятия при высокой биомассе } else {# Линейное увеличение от fmin до fmaxreturn(fmin + (fmax - fmin) * (b_ratio - blim) / (btarget - blim)) }}# Пример применения HCR для различных уровней биомассыb_ratios <-seq(0.2, 1.2, by =0.1)hcr_values <-sapply(b_ratios, harvest_control_rule)hcr_df <-data.frame(B_Bmsy = b_ratios, F_Fmsy = hcr_values)ggplot(hcr_df, aes(x = B_Bmsy, y = F_Fmsy)) +geom_line(linewidth =1.5, color ="blue") +geom_vline(xintercept =0.4, linetype ="dashed", color ="red") +geom_vline(xintercept =0.8, linetype ="dashed", color ="green") +labs(title ="Правило управления Harvest Control Rule (HCR)",x ="Отношение биомассы к Bmsy (B/Bmsy)",y ="Отношение fishing mortality к Fmsy (F/Fmsy)") +theme_minimal() +theme(plot.title =element_text(hjust =0.5))
# ------------------------- 7. МОДЕЛИРОВАНИЕ УПРАВЛЕНИЯ (MSE) ---------------------------# Упрощенное моделирование управленияmanagement_simulation <-function(init_biomass, r, k, catch, n_years =20, hcr_function) { biomass <-numeric(n_years +1) biomass[1] <- init_biomass catch_rec <-numeric(n_years) f_rec <-numeric(n_years) b_ratio_rec <-numeric(n_years)for (t in1:n_years) {# Расчет текущего отношения B/Bmsy b_ratio <- biomass[t] / (k /2) b_ratio_rec[t] <- b_ratio# Применение HCR для определения уровня изъятия f_multiplier <-hcr_function(b_ratio) f_rec[t] <- f_multiplier * (r /2) # F = multiplier * Fmsy# Расчет вылова на основе F catch_rec[t] <- f_rec[t] * biomass[t]# Обновление биомассы (модель Шефера) biomass[t +1] <-max(0.001* k, biomass[t] + r * biomass[t] * (1- biomass[t] / k) - catch_rec[t]) }list(biomass = biomass, catch = catch_rec, f = f_rec, b_ratio = b_ratio_rec)}# Запуск моделирования с различными начальными условиямиset.seed(123)n_sim <-10init_depletion <-runif(n_sim, 0.3, 0.7) # Различные начальные уровни деплетированностиsim_results <-list()for (i in1:n_sim) { init_biomass <- init_depletion[i] *median(msy_results$k) sim_results[[i]] <-management_simulation(init_biomass, median(msy_results$r), median(msy_results$k), Catch,hcr_function = harvest_control_rule)}# ------------------------- 8. АНАЛИЗ РЕЗУЛЬТАТОВ MSE ---------------------------# Подготовка данных для визуализацииyears_proj <-1:21sim_biomass <-sapply(sim_results, function(x) x$biomass)sim_catch <-sapply(sim_results, function(x) x$catch)# Визуализация траекторий биомассыbiomass_df <-data.frame(Year =rep(years_proj, n_sim),Biomass =as.vector(sim_biomass),Simulation =rep(1:n_sim, each =length(years_proj)))ggplot(biomass_df, aes(x = Year, y = Biomass, group = Simulation, color = Simulation)) +geom_line(linewidth =0.7) +geom_hline(yintercept =median(msy_results$b_msy), linetype ="dashed", color ="red") +labs(title ="Траектории биомассы в имитационном моделировании",y ="Биомасса") +theme_minimal() +theme(plot.title =element_text(hjust =0.5))
# Визуализация траекторий выловаcatch_df_proj <-data.frame(Year =rep(1:20, n_sim),Catch =as.vector(sim_catch),Simulation =rep(1:n_sim, each =20))ggplot(catch_df_proj, aes(x = Year, y = Catch, group = Simulation, color = Simulation)) +geom_line(linewidth =0.7) +geom_hline(yintercept =median(msy_results$msy), linetype ="dashed", color ="red") +labs(title ="Траектории вылова в имитационном моделировании",y ="Вылов") +theme_minimal() +theme(plot.title =element_text(hjust =0.5))
# ------------------------- 9. ФОРМУЛИРОВАНИЕ РЕКОМЕНДАЦИЙ ---------------------------# Формулируем рекомендации на основе анализаmedian_msy <-median(msy_results$msy)current_catch <-mean(Catch[(length(Catch)-2):length(Catch)])cat("АНАЛИЗ РЕЗУЛЬТАТОВ И РЕКОМЕНДАЦИИ:\n")
if (current_catch / median_msy >1) {cat("ВЫВОД: Текущий уровень вылова превышает оценку MSY. Рекомендуется сокращение вылова.\n")} elseif (current_catch / median_msy >0.8) {cat("ВЫВОД: Текущий уровень вылова близок к MSY. Рекомендуется осторожный подход.\n")} else {cat("ВЫВОД: Текущий уровень вылова ниже MSY. Возможно увеличение вылова.\n")}
ВЫВОД: Текущий уровень вылова близок к MSY. Рекомендуется осторожный подход.
cat("\nРЕКОМЕНДАЦИИ ПО УПРАВЛЕНИЮ ЗАПАСОМ:\n")
РЕКОМЕНДАЦИИ ПО УПРАВЛЕНИЮ ЗАПАСОМ:
cat("- Установить целевой уровень вылова в диапазоне ", round(median_msy *0.8, 2), "-", round(median_msy, 2), "тыс. тонн\n")
- Установить целевой уровень вылова в диапазоне 10.47 - 13.08 тыс. тонн
cat("- Внедрить правило управления (HCR) для адаптивного регулирования вылова\n")
- Внедрить правило управления (HCR) для адаптивного регулирования вылова
cat("- Усилить мониторинг запаса для улучшения оценок\n")
- Усилить мониторинг запаса для улучшения оценок
cat("- Проводить регулярную оценку запаса с обновлением рекомендаций\n")
- Проводить регулярную оценку запаса с обновлением рекомендаций
# =============================================================================# КОНЕЦ СКРИПТА# =============================================================================
17.1 Анализ результатов применения метода Catch-MSY: между надеждой и статистической реальностью
Что ж, давайте посмотрим, что у нас получилось в результате выполнения этого скрипта. Картина вырисовывается интересная, хотя и не без типичных для DLM сюрпризов.
Начнем с того, что наш метод Catch-MSY отработал как швейцарские часы — если, конечно, считать нормальным тот факт, что из 100000 итераций у нас осталось всего несколько тысяч жизнеспособных комбинаций параметров. Но это как раз и есть философия метода: мы не пытаемся найти единственно верное решение, а отбираем все биологически возможные сценарии.
Полученные оценки параметров говорят сами за себя. Коэффициент роста r со средним значением 0.14 — это вполне типично для многих промысловых видов. Хотя, если быть до конца честным, разброс значений от 0.01 до 1.5 и стандартное отклонение 0.13 намекают, что наша уверенность в этом параметре несколько преувеличена.
Емкость среды K оказалась в районе 555 тысяч тонн в среднем, что при медианном MSY в 13 тысяч тонн дает нам примерное представление о потенциале запаса. Любопытно, что медианное значение MSY практически совпадает со средним — 13.08 против 12.75 тысяч тонн, что намекает на относительно симметричное распределение.
А теперь самое интересное — текущий статус запаса. Средний вылов за последние три года составляет 12 тысяч тонн против оцененного MSY в 13.08. Отношение 0.92 — это тот самый момент, когда начинается настоящая игра в рулетку управления. Формально мы еще не превысили MSY, но мы уже в опасной близости от границы.
Но настоящий сюрприз ждал нас в оценке “деплетированности”. Медианное значение 3.13 — это, конечно, статистический нонсенс. Биомасса не может превышать емкость среды, если только мы не открыли новый закон популяционной динамики. Этот артефакт — классический пример ограничений метода Catch-MSY, когда при определенных комбинациях параметров модель выдает биологически невозможные значения.
Ориентиры управления выглядят более-менее разумно: MSY 13.08 тыс. тонн, Bmsy 243.8 тыс. тонн, Fmsy 0.048. Последнее значение особенно интересно — оно предполагает, что оптимальный уровень изъятия составляет около 5% от биомассы.
Имитационное моделирование управления показало, что при использовании предложенного правлила регулирования промысла (ПРП) средний вылов колеблется от 8 до 14 тысяч тонн в зависимости от начальных условий. Что характерно, даже в наихудших сценариях запас не коллапсирует — система управления оказывается достаточно устойчивой.
Графики траекторий биомассы и вылова — это настоящая поэма в данных. Мы видим, как различные сценарии сходятся к устойчивому состоянию около Bmsy, что собственно и является целью управления. Хотя некоторые симуляции показывают излишне консервативные подходы с выловом ниже потенциально возможного.
17.2 Заключение: между статистикой и управленческой реальностью
По итогам анализа можно сказать, что мы получили ровно то, что ожидали от метода Catch-MSY — ориентировочные оценки с изрядной долей неопределенности. MSY в районе 13 тысяч тонн выглядит разумной оценкой, учитывая историю промысла.
Текущий уровень вылова в 12 тысяч тонн — это тот самый момент, когда управляющим органам пора начинать тревожно задумываться. Мы еще не превысили лимит, но мы уже близки к той грани, где любое дополнительное давление может привести к перелову.
Рекомендация установить OДУ в диапазоне 10.5-13 тысяч тонн — это не проявление излишней осторожности, а единственно разумный подход в условиях высокой неопределенности.
Метод Catch-MSY в очередной раз доказал свою полезность как инструмент “первого приближения”. Он не дает нам точных оценок, но предоставляет достаточно информации для принятия взвешенных управленческих решений.
Ирония ситуации заключается в том, что чем меньше у нас данных, тем более сложные методы мы вынуждены использовать для их интерпретации. Но как показывает практика, иногда простые решения, основанные на ограниченной информации, оказываются более эффективными, чем сложные модели, требующие тонны данных.
В конечном счете, управление рыболовством — это не о поиске идеальных решений, а о выборе наименее плохих из доступных вариантов. И в этом смысле Catch-MSY предоставляет нам именно то, что нужно — достаточно информации, чтобы сделать этот выбор осознанно.
