ВведениеЧерноморско-азовская проходная сельдь Alosa immaculata (Bennett, 1835) –представитель полупроходных рыб Азово-Черноморского бассейна, где выделяют 2 популяции сельди, в зависимости от мест ее нереста, – в р. Дон (Донская популяция) и в р. Дунай (Дунайская популяция) [4]. Данное исследование выполнено для Донской популяции. Основной промысел донской сельди осуществляется в двух промысловых районах: в Керченском проливе – на путях нерестовой и зимовальной миграций через Азовское море в р. Дон (февраль-март и октябрь-ноябрь) ставными неводами и другими прибрежными сетными орудиями лова; и с апреля по июнь в р. Дон – закидными неводами [2; 20].Численность и биомасса сельди зависит от ряда абиотических и биотических факторов среды. Основными определяющими факторами являются: условия зимовки в Черном море; количество производителей, достигших мест нереста; величина весеннего разлива р. Дон; скорости течения, которые обеспечивают выживаемость икры сельди; температура и соленость воды; промысловая нагрузка. [11; 20]. Промысел сельди до конца 1990-х гг. осуществлялся целевыми орудиями лова – ставными сельдевыми неводами. После запрета на вылов сельди в 1998 г. и возобновление промысла в 2004 г., уловы сельди стали носить характер приловов в прибрежные орудия лова – ставные и закидные невода, вентеря, сетные орудия лова [15].Исторически оценка запаса сельди выполнялась при помощи метода прямого учета [1], по результатам учетной траловой съемки в Азовском море в летний и осенний периоды [7]. В настоящее время (начиная с 2019 г.) оценка запаса выполняется при помощи аналитических методов, рекомендованных институтом ФГБНУ «ВНИРО» [3]. В период до зарегулирования стока р. Дон (1930-1950 гг.) биомасса сельди составляла от 9,54 до 19,2 тыс. т [7; 15].После строительства и введения в эксплуатацию Цимлянского гидроузла в 1952 г., произошло перераспределение годового стока в р. Дон: сократился объем водного стока в весеннее время, снизились скорости течения, что привело к массовой гибели икры сельди, которая опускалась в придонные слои и погибала, не давая значительного пополнения. В результате, в указанный период, наблюдается резкое сокращение запаса сельди до 2,1-5,6 тыс. т [12].В 1970-1980 гг., в связи со строительством каскада низконапорных гидроузлов, запас сельди продолжил сокращаться. Биомасса в среднем за период составила 0,46 тонн. [2]. В 1994 г. запас сельди достиг абсолютного минимума и составил 0,15 тонн. В этот же год был введен запрет на промышленный лов сельди вплоть до 1996 г. [6; 8]. В последующий за запретным периодом год запасу сельди не удалось восстановиться и запрет специализированного промысла был продлен на период 1998-2004 гг. [7].  После возобновления промысла в 2004 г., в период 2007-2014 гг. наблюдалось постепенное увеличение объема общего годового вылова и умеренное восстановление запаса до уровня 3,2 тыс. т [7]. Стоит отметить, что сельдь – один из немногих видов, восстановивших свои запасы без искусственного воспроизводства [6; 15].В период 2014-2018 гг. промысловый запас сельди умеренно уменьшился за счет ряда факторов: низкий весенний водный сток, увеличение солености Азовского моря и Таганрогского залива, сокращение нагульных площадей молоди. В результате, в эти годы (кроме 2018 г. [7]) наблюдались неурожайные пополнения сельди. По результатам оценки прямыми методами учета, промысловая биомасса сельди в 2021 г. составляла 2,5 тонн.Основная цель работы – выполнение оценки промыслового запаса, промысловой смертности и целевых ориентиров эксплуатации популяции черноморско-азовской проходной сельди, по имеющимся данным в современный период в условиях недостаточности информационного обеспечения, при помощи методов математического моделирования.  Материалы и методыДля выполнения популяционного моделирования оценки биомассы запаса, промысловой смертности и биологических ориентиров использовались методы продукционного моделирования [26; 31], реализованные в современных программно-аппаратных комплексах в среде R [30]. Выбор в пользу ограниченных (упрощенных) продукционных моделей был обусловлен полнотой и обилием информационного обеспечения о популяции черноморско-азовской сельди. Когортные методы моделирования оказались не применимы по причине отсутствия непрерывных рядов данных из промысловых орудий лова – возрастной структуры, средней навески по возрастам из уловов и других данных. Применение полноценных продукционных моделей оказалось затруднительным по причине сложности стандартизации параметров уловистости различных орудий лова и индексов улова на единицу промыслового усилия. Высокая степень неопределенности ННН-добычи сделала затруднительным применение любых аналитических и статистических методов, основанных на информации о величинах промысловых уловов. Для решения проблемы высокой неопределенности в величинах суммарного годового промыслового изъятия применена методика косвенной оценки ННН-добычи [16; 17; 29; 33], ранее опубликованная авторами этой работы [18]. Выполнена оценка корректировочного коэффициента ННН-промысла – k. Данный коэффициент является поправочным коэффициентом к величинам официального вылова, позволяет косвенно учесть ННН-добычу, как долю к легальному вылову:ky=1+ICyCyгде: – ky — корректировочный коэффициент на долю ННН-промысла в год y;– ICy — незаконный улов в год y;– Cy — легальный улов в год y.Аппроксимация корректировочных коэффициентов выполнена за период 2000-2020 гг. с учетом ряда априорных допущений. В условиях полного запрета промысла и запрета специализированных орудий лова, для промысла сельди в 2000 г. [7; 15] было принято допущение об отсутствии специализированного браконьерского промысла – k2000 = 1. По результатам работы [19] было установлено, что уровень ННН-добычи в 2011 г. был не менее, чем в 2,7 раза выше, чем официальная добыча k2011 = 2,7. По результатам научного мониторинга структуры промысловых и научных уловов, в 2019 г. было отмечено, что незаконная добыча сельди, в качестве приловов, не менее чем в 3 раза превышает официальную статистику – k2019 = 3. Имеющиеся априорные оценки k2000, k2011, k2019 были аппроксимированы нелинейной логарифмической функцией для периода 2000-2020 гг., путем общего решения генерального дифференциального уравнения зависимости k от вектора времени t:dkdt=fkt=c+b*lntгде: - dk/dt – изменение корректировочного коэффициента k в непрерывном векторе времени dt;- fk(t) – интегрированная функция, описывающая изменение коэффициента k в точке во времени t;- t – момент времени;- c, b – параметры уравнения.Аппроксимация нелинейной функции и поиск параметров уравнения выполнен при помощи нелинейного моделирования в пакетах nls, nlstools в среде R [21]. Полученные оценки корректировочных коэффициентов, за период 2000-2020 гг., были переведены на шкалу уловов и сглажены методом полиномиальной регрессии для устранения межгодовой неопределенности – span = 0,3 [22]. Построение продукционной модели выполнено в два этапа. На первом этапе выполнялось построение упрощенной продукционной модели Шефера в реализации пакета CMSY [24; 25; 27], который ранее был доработан и апробирован для Азово-Черноморских видов [14]. При параметризации модели CMSY использовался ряд уточненной (на долю ННН-добычи) рыбохозяйственной статистики за период 2004-2020 годов. Параметр популяционной пластичности был определен как умеренный по таблице популяционной пластичности [28]. Период 2000-2003 гг. не использовался для моделирования по причине непромыслового статуса популяции (запрет промысла в эти годы) и крайне низких величин официального вылова, который был обусловлен по большей степени величиной научного изъятия. По результатам предварительного моделирования, при помощи CMSY, была построена модель при помощи комплекса JABBA [32]. Фреймворк JABBA позволяет выполнить расширенные процедуры стохастической диагностики модели и рассмотреть различные прогнозные сценарии, что качественно превосходит возможности модели CMSY. При параметризации диапазонов оптимумов в фреймворке JABBA использовались диапазоны оптимумов, полученные в результате оценки параметров модели CMSY (r, K, B0). Построение модели JABBA выполнялось в реализации «catchonly» – использовался только ряд скорректированных общих годовых уловов, без использования информации об улове на единицу промыслового усилия. Выполнялась апостериорная оценка ошибки процесса моделирования и расчет доверительных интервалов всех оцениваемых величин при p = 0,95. Прогнозирование по результатам моделирования выполнено с краткосрочным горизонтом 3 года, рассмотрены 6 сценариев, в зависимости от предполагаемого уровня промыслового изъятия: от 350 до 600 т с шагом 50 т для каждого сценария. РезультатыРезультаты расчета корректирующих коэффициентов на долю ННН-добычи, величины скорректированного общего годового улова представлены в таблице 1.  Таблица 1. Многолетняя рыбопромысловая статистика вылова черноморско-азовской сельди в водах России и УкраиныTable 1. Long-term fishing statistics of catch of Black Sea-Azov Pontic shad in the waters of Russia and UkraineГодОфициальный вылов, тАппроксимированные kОфиц. вылов + ННН сглаженный LOESS, т20000,2*10,0520010,752*1,471,3720023,102*1,746,1420037,781*1,948,42200410,82,0940,47200548,12,2191,66200648,82,31109,45200732,52,480,09200833,112,48101,9200970,422,55154,14201065,252,62187,92201169,842,68163,5201250,712,73170,3201385,732,78227,612014101,732,83249,6201578,542,87297,12016165,412,91411,162017189,172,95643,632018289,852,99762,882019233,5523,02776,892020230,443,05682,83* - официальный запрет на промысел, изъятие на научные нужды Рыбопромысловая статистика за период 2000-2003 гг. не использовалась при моделировании по причине крайне низких уловов и непромыслового статуса популяции. На основе ряда данных за 2004-2020 гг. (табл. 1), была построена предварительная модель CMSY [8]. Вычисленные параметры продукционной модели Шефера (табл. 2) использовались для априорной параметризации модели JABBA. Таблица 2. Результаты оценки параметров продукционного уравнения и ориентиров при моделировании на CMSYTable 2. The results of the evaluation of the parameters of the production equation and benchmarks in CMSY modelingПараметрОптимумДоверительный интервал при p = 0,95r0,5660,407 — 0,785K3432,91889,8B2005/BMSY1,05- На основе предварительной оценки параметров уравнения Шефера (табл. 2) выполнено её моделирование в фреймворке JABBA. Дополнительно при априорной параметризации использовались: - оценка ошибки в величине суммарного годового улова, catch.error = TRUE;- вариативность ошибки годового улова, catch.cv = 0,2;- оценка ошибки процесса моделирования, sigma.est = TRUE.Результаты оценки биомассы запаса (B), промысловой смертности (F) за период 2004-2020 гг. на фреймворке JABBA представлены в таблице 3. Относительные оценки биомассы запаса и промысловой смертности, к целевым ориентирам эксплуатации (BMSY, FMSY), представлены на рисунке 1. Были получены оценки целевых ориентиров эксплуатации: BMSY = 1975,8 т, FMSY = 0,29, MSY = 557,3 т и уточненные параметры r/K пары продукционного уравнения: r = 0,589, K = 3951,6 тонн. Таблица 3. Оценки биомассы запаса, промысловой смертности и их доверительных интервалов моделью JABBA, черноморско-азовская сельдь, 2004-2020 годыTable 3. Estimates of stock biomass, fishing mortality and their confidence intervals by the JABBA model, Black Sea-Azov shad, 2004-2020ГодBДоверительный интервал B,p = 0,95FДоверительный интервал F,p = 0,9520041315,5345,5 - 2755,40,060,01 - 0,1220051887,4686 - 3252,80,050,03 - 0,1320062479,9799,1 - 3670,50,040,03 - 0,1420072857,8952 - 3921,90,030,02 - 0,0820083173,31153,9 - 4201,40,030,02 - 0,0920093347,51358,7 - 4466,90,050,03 - 0,1120103395,51507 - 4702,60,060,04 - 0,1220113384,11612,2 - 4767,60,050,03 - 0,120123407,61702,1 - 4891,20,050,03 - 0,1201334111774,6 - 4925,50,070,05 - 0,1320143353,11783,8 - 4855,20,070,05 - 0,14201532941761,3 - 4835,30,090,06 - 0,1720163203,91718,5 - 4831,70,130,09 - 0,2420173012,41602,4 - 4645,80,210,14 - 0,4201826161276,8 - 4254,70,290,18 - 0,620192179,5810,9 - 3903,60,360,2 - 0,9620201791,3238,5 - 3667,60,380,19 - 2,86  Рисунок 1. Многолетняя динамика биомассы запаса относительно целевого ориентира (B/BMSY, рис. a) и промысловой смертности относительно целевого ориентира (F/FMSY, рис. b) и доверительные интервалы выполненных оценок (серая область) черноморско-азовской сельди по результатам моделирования на JABBAFigure 1. Long-term dynamics of stock biomass relative to the target (B/BMSY, Fig. a) and fishing mortality relative to the target (F/FMSY, (fig. b) and confidence intervals of the performed estimates (gray area) of the Black Sea-Azov herring based on the results of JABBA model ––В соответствии с полученными многолетними оценками и величинами целевых ориентиров эксплуатации, была построена многозональная схема правил регулирования промысла [23] (далее – ПРП), рисунок 2. Данная схема позволяет наглядно интерпретировать состояние популяции и влияние на нее промысловой смертности в различные годы.Рисунок 2. Многозональная схема правил регулирования промысла, продукционная кривая и целевые ориентиры состояния популяции и эксплуатации черноморско-азовской сельди по результатам моделирования. Красная область – коллапс популяции по причинам переэксплуатации; желтая область – буферная зона сокращения биомассы популяции без признаков переэксплуатации; оранжевая область – промысловая переэксплуатация популяции без признаков коллапса; зеленая зона – благоприятное состояние популяции без признаков переэксплуатацииFigure 2. A multi-zone scheme of fishing regulation rules, a production curve and targets for the state of the population and exploitation of the Black Sea-Azov shad based on the results of modeling. Red area – population collapse due to overexploitation; yellow area – buffer zone of population biomass reduction without signs of overexploitation; orange area – commercial overexploitation of the population without signs of collapse; green area – favorable state of the population without signs of overexploitation По результатам подгонки модели и на основе полученных оценок биологических и промысловых ориентиров, правил регулирования промысла, был выполнен краткосрочный прогноз динамики биомассы запаса и промысловой смертности, на период 2021-2023 гг., по различным сценариям промысловой нагрузки (рис. 3). Краткосрочный прогноз состояния биомассы запаса, при наиболее приемлемом рациональном сценарии, близком к MSY, представлен в таблице 4. Рисунок 3. Краткосрочный прогноз динамики биомассы запаса, относительно целевого ориентира (B/BMSY, рис. a), и промысловой смертности, относительно целевого ориентира (F/FMSY, рис. b), и доверительные интервалы выполненных оценок (серая область) черноморско-азовской сельди при различных сценариях промыслового изъятия в период 2021-2023 годыFigure 3. Short-term forecast of the stock biomass dynamics, relative to the reference point (B/BMSY, Fig. a), and fishing mortality, relative to the target reference point (F/FMSY, fig. b), and confidence intervals of the performed estimates (gray area) of the Black Sea-Azov shad under various scenarios of commercial fisheries in the period 2021-2023 Таблица 4. Сценарий краткосрочного прогноза биомассы запаса и промысловой смертнsости при суммарном годовом изъятии на уровне 450 т в период 2021-2023 годы Table 4. Short-term forecast scenario of stock biomass and fishing mortality with total annual withdrawal at the levl of 450 tons in the period 2021-2023ГодБиомасса запаса (B), тПромысловая смертность (F), доли202116760,285202217480,272202318210,259 ОбсуждениеРезультаты продукционного моделирования, полученные на JABBA, указывают на наличие высокой доли неопределенности и умеренной надежности полученных оценок. Высокая неопределенность во входных данных, вызванная неточностью косвенной оценки объемов ННН-добычи, существенно ухудшила надежность оценок модели, что в конечном итоге привело к широким доверительным интервалам полученных оценок. Тем не менее, несмотря на высокую неопределенность, полученные результаты позволили удовлетворительно описать современную многолетнюю динамику изменения биомассы запаса, промысловой смертности (табл. 3) и оценить ориентиры эксплуатации популяции сельди (рис. 1, 2). По результатам моделирования выявлено восстановление запаса сельди в период с 2004 по 2014 гг., что подтверждается результатами учетно-траловых съемок, проводимых в эти годы. В указанный период наблюдаются низкие показатели промысловой смертности, колеблющиеся в пределах от 0,03 до 0,07 (табл. 3), что дало возможность популяции сельди самостоятельно восстановить свои запасы без искусственного воспроизводства. В период 2014-2020 гг. отмечен постепенный тренд снижения биомассы запаса сельди, на фоне постоянно нарастающего тренда промысловой смертности, относительно предыдущего периода (с 0,07 в 2014 г. до 0,38 в 2020 г.). Начиная с 2018 г., уровень промысловой смертности стал превышать целевой ориентир FMSY = 0,28, что свидетельствует о переэксплуатации запаса сельди в период 2018-2020 годов.Детальное изучение состояния популяции, относительно целевых ориентиров (рис. 1,2), говорит о том, что практически за весь рассматриваемый период биомасса запаса сельди находилась выше целевого ориентира без признаков переэксплуатации. В период 2005-2017 гг. биомасса запаса находилась выше целевого ориентира в биологически безопасной зоне (B/BMSY >1), а промысловая смертность не превышала целевого ориентира эксплуатации (F/FMSY <1). Начиная с 2018 г., отмечены первые признаки переэксплуатации (F2018/FMSY ≈ 1), которая постепенно нарастает (F2020/FMSY = 1,3). Рост промысловой смертности в этот период отмечается на фоне постепенно снижающихся величин общего годового вылова: с учетом ННН – с 762 до 682 т., без ННН – с 289 до 230 т (табл. 1).Переэксплуатация в период 2018-2020 гг., на фоне сокращения величины общего годового изъятия (как с учетом ННН, так и без), обуславливалась постепенным трендом сокращения биомассы промыслового запаса (рис. 1, 2). Визуализация периодов изменения биомассы запаса и промысловой смертности графическим способом позволяет более детально изучить причинно-следственные связи (рис. 4).   Рисунок 4. Периоды трендов биомассы промыслового запаса (a) и промысловой убыли (b) черноморско-азовской проходной сельди по результатам моделирования. Точки – оценки биомассы запаса и промысловой убыли по результатам модели JABBA (табл. 3), линии – линейная регрессия оценок по вектору времени в выделенных интервалах, R2 – коэффициент детерминации линейной аппроксимации Figure 4. Periods of trends in biomass of commercial stock (a) and fishing mortality (b) of the Black Sea-Azov passing shad according to the results of modeling. The points are estimates of stock biomass and fishing mortality based on the results of the JABBA model (table 3), lines are the linear regression of estimates by the time vector in the selected intervals, R² is the coefficient of determination of the linear approximation Результаты трендового представления динамики биомассы запаса (рис. 4а) указывают на наличие 3-х различных периодов изменения биомассы: период восстановления биомассы запаса 2004-2009 гг., период стабильного, высокого уровня запаса 2010-2013 гг., период стремительного сокращения биомассы запаса 2014-2020 годы. Трендовая динамика промысловой смертности (рис. 4б) указывает на наличие двух периодов: период неопределенности промысловой смертности 2004-2013 гг., период стремительного роста промысловой смертности 2014-2020 годы. Анализ трендовых периодов совместно позволяет уверено утверждать, что сокращение биомассы запаса в период 2014-2020 гг.  напрямую обусловлено стремительным ростом промысловой смертности. Следует отметить, что в период 2004-2014 гг. промысловая смертность не являлась основным лимитирующим фактором биомассы запаса сельди.В настоящее время, 2019-2020 гг., популяция сельди находится на биологически небезопасном уровне с признаками незначительной, нарастающей переэксплуатации запаса (F2020/FMSY = 1,31). Вероятно, причиной переэксплуатации популяции после 2018 г. является постоянно нарастающий уровень ННН-добычи. Биомасса запаса, следующая монотонному линейному тренду на снижение с 2014 г. (рис. 4а), находится на уровне ниже целевого ориентира (B2020 = 1791,3 т., BMSY = 1975,8 т., B2020/BMSY = 0,9).Помимо увеличения промысловой нагрузки на запас сельди существенное влияние оказывают и абиотические факторы, в частности, сокращение водного стока в весенний период и затянувшееся осолонение Азовского моря [5; 13]. Зарегулирование стока р. Дон привело к значительному сокращению весеннего половодья, снижению скоростей течения реки, которые необходимы сельди для благоприятного нереста и, соответственно, поспособствовало увеличению солености в акватории моря [9; 10]. Изменившиеся условия окружающей среды стали причиной низкой эффективности размножения сельди, сокращения нагульного ареала молоди и ухудшения качественного и количественного состояния популяции в целом [10]. Учитывая отрицательный тренд биомассы запаса и состояния абиотической среды в 2014-2020 гг., в последующие годы, вероятнее всего, стоит ожидать дальнейшего сокращения объемов промыслового запаса сельди и ухудшение ее качественных характеристик.По результатам краткосрочного прогноза динамики биомассы запаса и промысловой смертности сельди на период 2021-2023 гг. был выбран сценарий, наиболее приемлемый с точки зрения рациональной эксплуатации по концепции MSY, на уровне добычи в 450 тонн. Этот сценарий был принят как базовый и представлен в таблице 4. Стоить отметить, что объем добычи в 450 т не является рекомендацией к установлению величины разрешенного вылова (РВ). Это теоретически возможный вылов с учетом ННН-нагрузки, доля которой, без принятия радикальных мер, составит около 300 т и не более 150 т останется на легальный вылов. Для более рациональной эксплуатации запаса нужно предпринимать радикальные меры борьбы с ННН-промыслом. До появления существенных сигналов сокращения объемов ННН-добычи рекомендуется сокращать объемы РВ до среднего уровня официальной добычи за последние 3 года – 250 тонн. При отсутствии мероприятий по предотвращению ННН-добычи, рекомендуется сохранять вылов на уровне 150-200 тонн. ВыводыРезультаты работы позволили получить оценки биомассы промыслового запаса, промысловой смертности и ориентиров эксплуатации за период исследования. Биомасса промыслового запаса находилась в пределах 1315-3411 т, промысловая смертность – 0,03-0,38. Целевые ориентиры состояния популяции: BMSY = 1975,8 т, FMSY = 0,29, MSY = 557,3 тонн.Выделены три периода состояния биомассы запаса: период восстановления 2004-2009 гг., период стабильного состояния 2010-2013 гг., период постепенного сокращения 2014-2020 годы. Выделены два периода изменения промысловой смертности: период неопределенности 2004-2013 гг., период постепенного увеличения 2014-2020 годы. Отмечено, что в период 2004-2013 гг. промысловая смертность не являлась ключевым лимитирующим фактором биомассы промыслового запаса сельди. В период 2018-2020 гг. постепенное увеличение уровня промысловой смертности привело к первым сигналам о незначительной переэксплуатации популяции. В 2020 г. переэксплуатация стала более выраженной, в результате чего биомасса промыслового запаса впервые опустилась ниже целевого ориентира эксплуатации. 5. Результаты работы обращают внимание на высокую долю ННН-промысла, который является существенной угрозой для популяции сельди. Авторы отмечают необходимость контроля и борьбы с незаконными добытчиками. 6. В настоящее время отсутствуют перспективы увеличения промысловой нагрузки на популяцию сельди. При осуществлении мероприятий по борьбе с незаконной добычей, допустим вылов на уровне 450 тонн. При отсутствии указанных мероприятий, рекомендуется сократить общий годовой вылов до 150-200 тонн.