CC BY
94
2008
Известия ТИНРО
Том 153
УДК 639.2.053.7
В.В. Суханов
Тихоокеанский научно-исследовательский рыбохозяйственный центр, Институт биологии моря ДВО РАН, г. Владивосток [email protected]
РАСЧЕТ ПРОМЫСЛОВОГО ЗАПАСА РЫБ ПО СОВОКУПНОСТИ ПРОБ
Предлагается метод расчета, позволяющий оценить промысловый запас популяции на акватории, покрытой множеством проб. Эта задача состоит из двух этапов: построение поверхности, интерполирующей пространственное распределение обилия в пробах, и численное интегрирование объема, накрываемого этой поверхностью. Дается описание двенадцати методов расчета таких поверхностей. Эти методы сравниваются с площадным методом Аксютиной. Обсуждается критерий, позволяющий выбрать метод, наилучший для данного материала. Излагаются основные операции с поверхностью: фильтрация, трансформация, бланкирование участков суши, расчет накрываемого объема, масштабирование. Метод иллюстрируется на примере популяции дальневосточной сардины в северо-западной части Японского моря.
Sukhanov V.V. Estimation of a fishery stock on the base of sample set // Izv. TINRO. — 2008. — Vol. 153. — P. 134-154.
New method is suggested allowing to estimate the fishery stock of population in the water area covered by a set of samples. The problem is solved in two steps: i) construction of the surface interpolating spatial distribution of the object abundance in the samples; ii) digital integration of the volume covered by this surface. Twelve different approaches for this surface construction are presented. The results of the stock estimations with these approaches are compared with the common Aksjutina "area method". The criterion is discussed for choice of the best method in dependence on data properties. Basic operations with the surface are described, as filtration, transformation, blanking of land areas, calculation of the covered volume, scaling. The method application is illustrated on the example of the stock estimation for japanese sardine in the northwestern Japan Sea.
Введение
Задача расчета запаса нередко возникает в рыбопромысловых исследованиях. Ее можно сформулировать следующим образом. На акватории осуществляется сбор стандартизированных проб. Каждая проба представляет обилие некоего вида рыб в данной точке сбора. По этим материалам требуется определить общий промысловый запас, т.е. суммарное по всей акватории обилие той части популяции, размер тела особей которой не ниже установленного уровня.
Эта задача в самой элементарной форме была сформулирована планктонолога-ми Гензеном и Апштейном (Hensen, Apstein, 1897, цит. по: Никольский, 1965). Потом она неоднократно видоизменялась и уточнялась. В окончательном виде способ решения этой задачи получил название "площадного метода" (Аксютина, 1968).
Площадной метод можно выразить в виде формулы: Z = zS / (ks), где Z — суммарный промысловый запас популяции, z — средний улов в одной пробе, S — площадь всей акватории, s — средняя площадь зоны облова для одной пробы, k — коэффициент уловистости орудия лова. В уловы входят только те особи, размер тела которых оказался не менее установленного промыслового размера. Статистическая ошибка площадного метода в наибольшей степени определяется ошибкой параметра z — среднего улова в одной пробе. Эта ошибка выражается отношением sj\iN, где s — среднеквадратическое отклонение уловов по всем N пробам. Кроме описанного базового варианта З.М. Аксютина (1968) ввела в свой метод модификации. Она предложила разделять акваторию на однородные участки (метод изолиний) и вместо средней арифметической использовать среднюю геометрическую для оценки величины z.
Основная специфика классического площадного метода состоит в том, что пробы считаются пространственно независимыми друг от друга и имеют одинаковый статистический вес. Эти особенности были расценены ихтиологами как недостатки метода. Поэтому ими были изобретены различные усовершенствования, которые обобщены в так называемом "методе страт" (Столяренко, Иванов, 1988; Борисовец и др., 2003), который фактически представляет собой развитие метода изолиний Аксютиной. Перед расчетами вся акватория разбивается на совокупность однородных участков, которые называются "стратами". Каждая страта характеризуется своим уровнем популяционного обилия и занимает конкретный ареал на акватории. Для оценивания среднего улова z в одной пробе вместо простой невзвешенной средней используется оценка, средневзвешенная по всем стратам. При этом веса таких страт пропорциональны занимаемым ими площадям на изучаемой акватории.
Цель настоящей работы заключается в дальнейшем развитии как классического площадного метода, так и его средневзвешенных модификаций. При этом главное внимание мы будем уделять компьютерным картографическим методам работы с цифровой моделью поверхности. Эта поверхность будет описывать пространственное распределение популяционного обилия. Часть статьи с описанием этих картографических методов представляет собой краткие выжимки текста из руководства В.В. Суханова (2005). На такое повторение пришлось пойти потому, что "Известия ТИНРО" имеют гораздо более широкий круг читателей, чем книжка, выпущенная в местном издательстве малым тиражом.
Все расчеты в этой работе были осуществлены с помощью программы Surfer, ver. 8, разработанной фирмой Golden Software Inc, USA. Именно на основе этой эффективной программы рекомендуется решать задачу расчета запаса. Однако это не обязательно, если сходные инструменты найдутся в других программах. В частности, аналогичные задачи можно решать при помощи программ, обеспечивающих эксплуатацию геоинформационных систем.
Результаты и их обсуждение
Методы расчета сетки
Вначале нужно построить цифровую модель пространственного распределения уловов. Пусть х — географическая долгота, y — географическая широта, z — величина улова в стандартной пробе с координатами x, y. Расчет цифровой модели поверхности, описывающей пространственное распределение обилия популяции, приводит к тому, что на базе точек, случайно распределенных по акватории, создается регулярная сетка, в ху-узлах которой располагаются z-значения уловов.
Характеристика "случайно распределенные" означает, что точки наблюдений могут быть беспорядочно, незакономерно распределены по плоскости карты. При этом расстояние между соседними точками наблюдений не обязано сохра-
135
няться постоянным. Вместе с тем цифровая модель поверхности z = f(x, y) требует, чтобы xy-координаты точек были регулярно размещены на плоскости, т.е. изменялись по х- и по y-координатам с постоянным шагом.
Когда xy-координаты наблюдений беспорядочно распределены по карте, между ними и за ними существуют незаполненные пустоты. При расчете регулярной сетки эти пустоты заполняются интер- или экстраполированными значениями z-координаты, которая обозначает величины уловов в стандартных пробах. Программа Surfer ver. 8 поддерживает 12 методов построения равномерной (регулярной) сетки. В каждом методе z-значение сеточного узла вычисляется по специфическому алгоритму и может приводить к своей количественной интерпретации исходных данных.
Большинство методов расчета регулярной сетки используют интерполяционные алгоритмы с взвешенной средней. Это значит, что при прочих равных условиях, чем ближе к узлу сетки расположена данная эмпирическая точка, тем большим становится ее статистический вес при расчете z-координаты этого узла. Обсудим особенности каждого метода. При ссылке на методы будут использоваться как русские, так и английские их названия (Суханов, 2005). Это облегчит пользователю ориентацию в списках меню программы Surfer.
Вид поверхности, порождаемой разными, сильно отличающимися друг от друга методами, будет проиллюстрирован на одном и том же исходном фрагменте поля. Эти поверхности намеренно будут рассчитаны стандартным способом, без особой настройки параметров, т.е. при тех значениях опций, которые приняты для используемых методов по умолчанию.
Inverse Distance to a Power. Метод "обратно пропорционально расстоянию в некоторой степени" — это очень быстрый метод расчета сетки. Вес, задаваемый конкретной точке данных при вычислении узла сетки, пропорционален величине, обратной расстоянию от точки до узла, причем это расстояние возведено в некоторую степень. При вычислении узла сетки эти веса пересчитываются в доли таким образом, что их сумма становится равной единице. Когда xy-координаты точки в точности совпадают с узлом сетки, статистический вес этой точки приравнивается к единице, а веса всех остальных точек становятся равными нулю. Иными словами, данный метод может работать как точный интерполятор.
Степенной параметр управляет тем, как именно уменьшается вес по мере роста расстояния от точки наблюдения до узла сетки. Чем больше эта степень, тем больший вес имеют близкие точки; чем меньше степень, тем равномернее распределяются веса среди точек. По умолчанию эта степень принята равной двойке.
Одно из досадных свойств данного метода состоит в том, что он порождает "бычьи глаза" (локальные круглые пятна на карте, окружающие точки наблюдений). На рис. 1 показан пример с фрагментом поверхности, рассчитанной при помощи обсуждаемого метода. Хорошо видны многочисленные ямки вдоль ложбин и холмики вдоль гребней. Каждый из этих артефактов расположен в точке,
где была взята отдельная проба. На карте изолиний эти ямки и холмики превратятся в светлые и темные "бычьи глаза".
Рис. 1. Фрагмент поверхности, рассчитанной при помощи метода "обратно пропорционально расстоянию в некоторой степени" Fig. 1. Surface fragment computed with the help of "Inverse D is-tance to a Power" method
Для того чтобы такие артефакты стали слабее или исчезли вообще, здесь можно установить ненулевое значение для параметра сглаживания (Smoothing). Такой параметр приводит к тому, что даже если узел сетки находится точно в точке наблюдения, этой точке не присваивается абсолютно вся единичная весовая доля. Параметр сглаживания сокращает эффект "бычьих глаз" путем как бы размазывания интерполированной сетки. В таком случае данный метод переходит в категорию сглаживающих интерполяторов.
Kriging. В целом геостатистика как совокупность математически обоснованных методов впервые появилась в геологии в связи с оценкой запасов руд. Отдельные методы применялись и ранее, в частности горным инженером Криге (Мальцев, 1999), в честь которого и назван основной расчетный метод геостатистики — кригинг. Впоследствии кригинг стал применяться и в других областях, где требуется аккуратный расчет полей.
Кригинг пытается учесть и сохранить в рассчитываемой сетке тенденции, которые содержатся в исходных данных: так, например, чтобы высокие точки, расположенные вдоль линии, могли быть связаны в гребень, а не выглядели цепочкой отдельных пятен типа "бычий глаз".
На рис. 2 показан фрагмент поверхности, рассчитанной при помощи данного метода. Отчетливо видны три параллельных гребня, разделенные двумя ложбинами. Вместе с тем на поверхности, хотя и гораздо слабее, чем на рис. 1, но все-таки также заметны артефакты в местах расположения эмпирических точек: мелкие проколы в ложбинах и невысокие пички на гребнях. Эти и другие артефакты можно убрать, если сделать специальную настройку параметров. В кригинг включены три группы настраиваемых параметров: тип дрифта, модель вариограммы и эффект самородка.
Рис. 2. Фрагмент поверхности, рассчитанной при помощи кригинга
Fig. 2. Surface fragment computed with the help of kriging
Тип дрифта (Drift Type) — это важный параметр в дополнительных опциях кригинга. Когда точки с данными равномерно рассеяны внутри области, выбор типа дрифта не влияет на рассчитываемую там сетку. Этот выбор становится важным, когда нужно интерполировать внутри протяженной пустой области или экстраполировать за пределы области, охваченной данными.
В программе Surfer доступны три варианта дрифта: None (без дрифта), Linear (линейный), Quadratic Drift (квадратичный дрифт). Когда есть сомнения, лучше всего использовать тип None, предполагая, что ваши данные рассеяны на плоскости более или менее равномерно. Если данные варьируют в соответствии с линейной, монотонной тенденцией (на одном краю акватории точки наблюдений распределены гуще, чем на противоположном краю), тогда целесообразно выбрать тип Linear Drift. Если же данные варьируют согласно квадратичной тенденции (например, параболическая чаша: в центре акватории точки рассыпаны гуще, чем на периферии, или наоборот — в центре реже, чем по краям), тогда лучше всего использовать Quadratic Drift.
Вариограмма (Variogram) есть количественная мера того, насколько быстро z-переменная изменяется в пространстве. Две точки с близкими xy-координа-тами в среднем будут иметь более похожие значения их z-координат, чем две
137
далеко отстоящие друг от друга точки. Вариограмма отражает степень разброса по z для двух произвольно выбранных точек поля по мере увеличения дистанции между ними. На поверхности, изрезанной крутыми перепадами уровней уловов, степень различий между двумя пробными точками растет быстрее, чем на плавной, медленно меняющейся поверхности. В реальных данных часто присутствует анизотропность: например, по мере удаления от берега биомасса прибрежного вида неуклонно снижается. Поэтому в вариограмме степень различий между пробными точками в общем случае зависит не только от дистанции между пробными точками, но и от направления в пространстве. Программа Surfer ver. 8 поддерживает 12 самых известных моделей вариограммы и позволяет оценить параметры каждой из них. Параметры вариограммы используются для расчета весов наблюдений при расчете текущего узла сетки.
Эффект самородков (Nugget Effect) также представляет собой важный параметр в кригинге. Его значение может изменяться от нуля до общей дисперсии по всему полю. В геостатистике эффектом самородков называют ненулевое численное значение вариограммы при нулевом расстоянии между точками. Большая чем нуль величина эффекта самородков говорит о неоднозначности поля: в некоторых местах это поле представлено несколькими несовпадающими значениями z-координаты при дублирующих ху-координатах. Само название этого термина указывает на неоднозначность поля: в данной точке обнаружилась пустая порода и содержание золота равно нулю, а буквально рядом попался самородок и это содержание подскочило до очень большой величины.
В распределениях уловов эффект самородков возникает в том случае, если исходные данные содержат ошибки. Если повторные измерения улова z в одной и той же точке на ху-плоскости осуществляются со случайными инструментальными погрешностями (из-за особенностей орудия лова) или находятся под влиянием естественных природных флюктуаций, то эффект самородков становится явно выраженным. Проявление эффекта самородков зависит от выбранной модели вариограммы. После его включения в кригинге усиливаются свойства сглаживающего интерполятора: больше доверия общей тенденции, а не индивидуальным точкам. Чем выше заданная величина эффекта самородков, тем глаже итоговая сетка.
Minimum Curvature. Метод минимальной кривизны широко используется в науках о Земле. Интерполированная поверхность, порожденная этим методом, аналогична тонкой линейно-эластичной плоскости, проходящей через каждую из эмпирических точек с минимально возможной величиной искривления. Примерно так в статистике работают так называемые кубические сплайны. Этот метод интерполяции имитирует технику чертежников начала 20-го века. Для того чтобы начертить плавную кривую, проходящую по цепочке точек, они втыкали в эти точки иголки, потом помещали между ними, изогнув, упругую узкую стальную полоску (именно она называлась сплайном) и использовали ее как лекало для черчения кривой.
Вместе с тем этот метод не относится к группе точных интерполяторов. Метод минимальной кривизны действует примерно так же, как если бы некто пытался накрыть эмпирические точки мозаикой, сшитой из маленьких пластинок упругого металла. Чем выше внутреннее напряжение в этих пластинках (параметр Internal Tension), тем сильнее они сопротивляются изгибу, тем больше аппроксимируемая поверхность напоминает черепицу из жестких фасеток. Параметр Boundary Tension также управляет величиной напряжения, но теперь уже не в серединных областях поля, а на краевых, граничных его участках.
Modified Shepard's Method. Усовершенствованный метод Шепарда использует метод Inverse Distance to a Power, взвешенный методом наименьших квадратов. При локальном использовании метода наименьших квадратов феномен "бычьих глаз" сильно ослабляется или вообще исчезает. Параметр Smoothing
138
Factor позволяет методу Шепарда работать не только как точному, но и как сглаживающему интерполятору. Ч ем больше значение этого параметра, тем сильнее эффект сглаживания. Наиболее разумно значение, лежащее между нулем и единицей.
Triangulation with Linear Interpolation. Триангуляция с линейной интерполяцией предельно бережно относится к исходным эмпирическим данным. Метод работает путем создания треугольников, стороны которых — это прямые отрезки, соединяющие точки взятия проб. Эти точки связываются таким образом, чтобы ни один из созданных треугольников не пересекался другими. В результате получается поверхность, склеенная из смежных треугольных фасеток. Каждый треугольник определяет плоскость, проходящую через его вершины, которые находятся непосредственно в точках с исходными данными. Поэтому этот метод является точным интерполятором.
Разбивать плоскость на треугольники можно разными способами даже для одних и тех же данных. Лучше всего поверхность поля аппроксимируется фасетками, как можно более близкими к равносторонним треугольникам. Такое разбиение плоскости дает оптимальная триангуляция Делоне.
На рис. 3 показан тот же самый фрагмент поверхности, что и на рис. 1, 2, но рассчитанный теперь при помощи триангуляции с линейной интерполяцией. Гребни и ложбины хорошо проявились, ямки и холмики исчезли, но вся поверхность приобрела какой-то изломанный, скомканный характер.
Рис. 3. Фрагмент поверхности, рассчитанной при помощи триангуляции с линейной интерполяцией
Fig. 3. Surface fragment computed with the help of triangulation and linear interpolation
Триангуляция работает лучше, когда имеется достаточное количество пробных точек, более или менее равномерно распределенных по акватории. Точки, далеко отстоящие от общей группы, порождают резко выделяющиеся, узкие треугольные фасетки. Если число точек менее 200, то фасетки также выглядят слишком выделяющимися.
Natural Neighbor. Метод естественного соседа базируется на сети полигонов Тиссена, которые, в свою очередь, строятся на основе классического метода триангуляции Делоне. Построение полигональных областей Тиссена приводит к следующему результату: каждая точка, обозначающая траловую пробу, окружается замкнутой ломаной границей — такой, что все точки внутри этой границы оказываются ближе к данной точке с траловой пробой, чем к каким-либо другим точкам акватории (Sibson, 1980; Ripley, 1981). Границы полигонов Тиссена являются отрезками перпендикуляров, восстановленных к серединам сторон треугольников в триангуляции Делоне (рис. 4).
Интерполяционный алгоритм Natural Neighbor при расчете очередного узла сетки использует взвешенную среднюю по соседним наблюдениям, в которой веса пропорциональны площадям этих соседних полигонов. Данный метод не экстраполирует контуры, выходящие за границы области, которая образована выпуклой оболочкой, натянутой на фактические точки. Иными словами, этот метод не выходит за границы полигонов Тиссена.
Рис. 4. Построение полигонов Тиссена (сплошные линии) по треугольникам Делоне (штриховые линии)
Fig. 4. Construction of Thiessen polygons (solid lines) on the base of triangles of Delauney (hatch lines)
Nearest Neighbor. Метод ближайшего соседа назначает для текущего узла сетки значение z-координаты, взятой из ближайшей к этому узлу эмпирической точки. Этот метод имеет много названий, одно из которых, "диаграмма Вороного", было использовано Е.Э. Борисовцом с соавторами (2003). Метод ближайшего соседа родственен упомянутому выше методу естественного соседа. Вместе с тем здесь при вычислениях узлов сетки не используется средняя взвешенная по площадям соседей, поэтому рассчитанная поверхность получается как бы ступенчатой (рис. 5). В центре торца каждого столбика находится эмпирическая точка, определяющая в своей ближайшей окрестности зафиксированный для нее
уровень уловов. Хотя это тот же самый фрагмент поверхности, что и на рис. 1, 2, 3, выглядит он не совсем похоже, скорее напоминая неаккуратно уложенную брусчатку.
Рис. 5. Фрагмент поверхности, рассчитанной при помощи метода ближайшего соседа
Fig. 5. Surface fragment computed with the help of "Nearest Neighbor" method
Этот метод эффективен для выявления пропусков в данных. Иногда встречаются почти полные и почти равномерные сетки. Но в них наблюдаются несколько пустых областей, которые требуется исключить из рассчитываемой сетки из-за отсутствия в этих областях реальной информации. Для этого нужно установить достаточно малые радиусы для эллипса поиска (Search Ellipse). Тогда пустые области в создаваемой сетке будут забланкированы (оставлены пустыми).
Если представить пространственное распределение уловов таким, как оно "видится" с точки зрения площадного метода Аксютиной, то это будет нечто, похожее на рис. 5. При этом все столбики будут иметь одинаковую площадь сечения и могут быть в произвольном порядке перетасованы по плоскости: для средней арифметической безразлично, где именно была взята та или иная проба, поскольку предполагается, что все они независимы друг от друга.
Polynomial Regression. Метод полиномиальной регрессии хорош, когда нужно выявить глобальные пространственные тенденции в фактических данных.
Он не является истинным интерполятором, потому что не пытается точно воспроизводить наблюдаемые z-значения уловов. Подгонка коэффициентов полинома осуществляется сразу по всем имеющимся в наличии эмпирическим точкам. Численные значения для коэффициентов подогнанного полинома выводятся в специальный Grid Report (отчет по расчету сетки).
Группа опций Surface Definition позволяет задать желаемый тип полиномиальной регрессии. Возможны пять вариантов, среди которых наиболее богат настройками последний вариант User defined polynomial — заданный пользователем полином. Группа опций Parameters позволяет задать максимальные степени для х- и y-переменных в уравнении полинома. Все установки тут же отображаются в форме записи полиномиального уравнения. Оценки коэффициентов этого уравнения можно прочитать в разделе Gridding Rulers в рапорте по итогам расчета сетки Gridding Report. Чтобы вывести этот рапорт на экран, надо включить флажок в окошке Grid Report в настройках метода.
Radial Basis Functions. Эта группа включает различные методы интерполяции данных, которые являются точными интерполяторами, поскольку строят поверхность, проходящую прямо по точкам. Вместе с тем здесь также можно ввести сглаживающий фактор, порождающий более пологую поверхность.
Радиальные базисные функции, которые нужно определить, аналогичны ва-риограмм-моделям в методе кригинга. Эти функции определяют оптимальное множество весов, "навешиваемых" на эмпирические точки в процессе их интерполяции в узел сетки. Среди возможных настроек отметим следующие важнейшие.
Опции из группы Basis Function позволяют определить параметры функций для операции расчета сетки. Список типов содержит перечень базисных функций, используемых при вычислении узлов сетки. Для большинства случаев лучше всего подходит мультиквадратическая функция (Multiquadric function).
Константа R_Parameter — это сглаживающий параметр. Чем больше эта величина, тем округлее получаются вершины и ямы поверхности и тем более гладкими становятся изолинии.
Moving Average. Метод скользящей средней вычисляет z-значение в текущем узле сетки путем простого усреднения всех точек, попавших в эллипс поиска, центр которого находится в этом узле. В разделе меню Advanced Options можно определить, каково должно быть минимально допустимое число точек, попадающих в этот эллипс. Если оно окажется меньшим, то данный сеточный узел будет забланкирован.
Data Metrics. Метод метрики данных собирает информацию о данных и создает регулярную сетку, описывающую эту информацию. Локальное подмножество данных, используемое для сбора такой информации, определяется параметрами эллипса поиска (Search). При расчете очередного узла привлекаются следующие локальные (в окрестности данного узла) статистические показатели: ранжированные z-статистики; моменты для координаты z (средняя, дисперсия и т.д.); прочие статистики по z; сведения по xy-локализациям точек; топографические данные.
Local Polynomial. Метод локального полинома вычисляет z-координату текущего узла сетки путем полиномиальной аппроксимации всех точек, попавших в эллипс поиска (Search). Можно задать геометрические параметры этого эллипса и минимально допустимую заполненность его точками. Там же можно установить величину (не очень большую) для порядка степени у локального полинома.
Перекрестная проверка. Команда меню Grid | Residuals вычисляет разность между каждым z-значением по фактическим данным и между соответствующим ему интерполированным z-значением на сеточной поверхности. Корень из среднего квадрата этих остатков есть количественная мера того, насколько хорошо данная сетка аппроксимирует реальные данные. Вместе с тем такой способ
не вполне чистоплотен: при расчете сетки используются те же данные, которые потом привлекаются для проверки. Честнее было бы привлекать для тестирования независимую информацию. Именно такой подход используется в методе перекрестной проверки.
Кнопка Cross Validate ... позволяет сделать настройку перекрестной проверки сетки на соответствие фактическим данным. Эта кнопка доступна в основном диалоговом окне Grid Data.. Окно открывается сразу после того, как, выбрав команду меню Grid | Data, указать файл с исходными данными для расчета сетки. Метод перекрестной проверки можно использовать для сравнения интерполяции xyz-данных при помощи разных методов расчета сетки. Итоги этой проверки можно также применять при огрубленном оценивании качества аппроксимации в разных участках карты — для планирования мест дополнительного сбора эмпирических данных.
Ошибка интерполяции, порождаемая выбранным методом расчета сетки, вычисляется как разность между интерполированным значением z и наблюдаемым значением. При этом перед интерполяцией в данную точку само наблюдаемое значение z-координаты временно изымается из таблицы с данными, чтобы не влиять на результат собственной проверки. Иными словами, для того чтобы протестировать тот или иной метод расчета сетки на материале из N эмпирических точек, нужно N раз рассчитать эту сетку, каждый раз вначале изымая из данных очередную эмпирическую точку. Это та самая точка, z-координата которой в настоящий момент подвергается тестирующей интерполяции при помощи сеточной поверхности, построенной по всем остальным N—1 точкам. Идея используемого здесь подхода в некоторой степени напоминает метод Кенуя, уменьшающий смещение оценок при подгонке параметров нелинейных моделей (Quenouille, 1956).
Операции с поверхностью
Поверхность, аппроксимируемая равномерной сеткой, может быть подвергнута дополнительным манипуляциям. С реди них пять важнейших операций: фильтрация шумов, математическое преобразование сетки, бланкирование участков суши, масштабирование и та операция, ради которой первоначально и задумывалась эта статья — вычисление суммарного объема, находящегося под рассчитанной поверхностью.
Фильтрация шумов. На поверхности обычно наблюдается высокочастотный шум: резкие и незакономерные скачки в величине рассчитанных уловов. Этот шум следует погасить при помощи процедуры фильтрования. Обычно фильтр представляет собой квадратную матрицу с нечетным числом строк и столбцов. Числа этой матрицы суть веса узлов, расположенных в маленьком квадратном фрагменте сетки. Веса убывают от центрального узла к краям матрицы. Эти веса участвуют при расчете улова, средневзвешенного по узлам сетки, которые накрываются матрицей фильтра. Полученная средневзвешенная оценка помещается в тот узел исходной сетки, который совмещен с центральным, срединным узлом матрицы. После такого усреднения величина улова в центральной точке фильтра теперь уже отчасти находится под влиянием уловов в соседних узлах матрицы. Если прежде в этой центральной точке наблюдался резкий скачок величины улова, то после фильтрования он частично погашается за счет влияния соседей. Затем центр матрицы фильтра перемещается на следующий узел исходной сетки, и операция взвешенного усреднения повторяется заново.
Отфильтрованную сетку можно фильтровать повторно и тем самым добиваться дальнейшего понижения в уровне высокочастотных шумов. Если исходная сетка рассчитана с мелким шагом (число строк и число столбцов достигает многих сотен), то высокочастотные шумы за один проход погашаются слабо и операцию фильтрования приходится повторять многократно. Не существует об-
щепринятых правил, сколько проходов фильтрации нужно делать и какой надо выбрать сглаживающий фильтр (а их насчитывается около 50 разновидностей). Можно предложить следующую рекомендацию.
Надо вычислить не сглаженную фильтром сетку, а потом посмотреть карту, полученную на ее основе. Обычно на такой карте видны узкие пики и проколы наподобие изображеных на рис. 1 или 2. Эти артефакты указывают на то, что при настройках метода расчета сетки не удалось установить достаточный уровень сглаживания аномальных выбросов, содержащихся в данных. Такие дефекты устраняются фильтрацией сетки. Для начала можно выбрать фильтр, предлагаемый в программе Surfer по умолчанию — это Gaussian (3x3). После просчета каждого варианта фильтрации нужно просматривать результаты его работы на карте. Число проходов фильтра (Number of Passes) сначала следует установить небольшим, а потом, если высокочастотные флюктуации с карты еще не исчезли, увеличивать значение этого параметра.
Математическое преобразование узлов. Выбирается каждый узел сетки, и его z-значение величины улова подвергается одному и тому же математическому преобразованию. Это преобразование может быть как совсем простым, так и достаточно сложным. Например, в простом случае нужно выразить уловы в других единицах измерения или сделать поправку на коэффициент уловисто-сти — для этого каждое z-значение сетки надо пересчитать с учетом масштабного коэффициента. В программе Surfer есть команда Grid | Math..., активизация которой открывает строку текстового редактора, где можно набрать формулу преобразования, которое нужно проделать с сеткой. Приведем несколько примеров.
Нередко после расчета сетки в некоторых ее узлах появляются отрицательные значения уловов (результат неадекватной интер- или экстраполяции). Эти абсурдные значения нужно заменить нулями, для чего в окне текстового редактора набрать формулу C = max(0,A), где исходная сетка обозначается буквой А, а итоговая сетка — буквой С. Математическая операция взятия максимума max(0,A) работает следующим образом. Если конкретный узел в сетке А имеет отрицательное z-значение улова, то это значение заменяется нулем. В противном случае, когда z > 0, оно оставляется неизменным. Попутно отметим, что метод Inverse distance to a power в принципе не может порождать z-значения, выходящие за минимальную и максимальную границы, обнаруженные в эмпирических данных. То же самое можно сказать и про метод Moving Average.
Возможны не только унарные преобразования сетки, но и бинарные операции с двумя сетками. Так, если требуется сложить z-значения из двух сеток А и В, а результат записать в сетку С, то в окне текстового редактора нужно набрать формулу С = А + В. Сетки А и В должны иметь одни и те же геометрические характеристики (габариты и шаг между узлами).
Можно создать новую сетку, которая будет рассчитана не на базе фактических измерений, а на основе математической формулы. Так, например, рассмотрим условно-практическую задачу. На исследуемой акватории нужно выделить перспективные участки, в пределах которых прибрежное рыболовство наиболее выгодно с экономической точки зрения. Наибольшая прибыль получается в том случае, когда разность между доходами и расходами достигает максимального значения. Эту задачу, если сформулировать ее в виде учебного примера, т.е. предельно просто и кратко, можно решить следующим образом.
Сначала посчитаем расходы. Предположим, что ежедневно после окончания лова судно возвращается к пункту приема сдавать улов. Будем считать, что себестоимость добытого и доставленного на берег улова слагается из постоянной части, равной 3000 руб./ сут и переменной части. Эта переменная часть расходов линейно зависит от расстояния между местом промысла и пунктом приема улова и составляет 100 руб. на один километр пути.
143
Точка, где расположен пункт приема рыбы, представлена x-координатой, равной 135 км, а также y-координатой, равной 45 км. В этом условном примере предполагается, что начало координат (x = 0 км, y = 0 км) расположено в левом нижнем углу картосхемы, изображающей акваторию района промысла. Тогда расстояние между пунктом приема рыбы и местом ее промысла с координатами
(x, y) выражается функцией x-135) + (y-45) . Рассчитаем сетку, описывающую пространственное распределение суточных расходов на обеспечение вылова рыбы и на доставку улова к месту сдачи. Эти расходы описываются суммой,
состоящей из постоянной и переменной частей: 3000 +100 -^(x -135)2 + (y - 45). Активизируем меню Grid | Function, и в окне текстового редактора наберем формулу Z = 3000 + 100*sqrt(pow((X-135),2) + pow((Y-45),2)). Оператор sqrt(...) представляет функцию квадратного корня, а оператор pow(...,2) обозначает функцию возведения величины ... в степень, равную двум. В итоге получаем сетку (в дальнейшем она будет обозначена буквой B), представляющую пространственное распределение расходов, т.е. распределение себестоимости улова.
Перейдем теперь к доходам. Рассчитаем сетку A, описывающую распределение прибрежных уловов по акватории. Уловы будем выражать в килограммах рыбы, которая вылавливается в течение рабочего дня. Допустим, что в пункте приема цена одного килограмма пойманной рыбы составляет 15 руб. Если пересчитать сетку уловов A с помощью выражения 15*A, то мы получим сетку, которая описывает пространственное распределение потенциальных доходов от рыбы, вылавливаемой за день.
Теперь можно построить сетку, описывающую пространственное распределение прибыли от вылова рыбы. Командой Grid | Math. вызываем окно текстового редактора, куда впечатываем формулу C = 15*A—B, где сетка А содержит уловы, сетка 15*A — доходы, сетка B — расходы, а сетка С — прибыль от вылова и продажи рыбы (налогами пока пренебрегаем). Сделаем последнюю операцию: нарисуем карту по сетке С. Те места на акватории, где ожидаются высокие уровни прибыли, следует считать самыми выгодными сейчас для промысла.
Перечень встроенных функций, которые знает программа Surfer, довольно длинен. Он представлен в меню Help | Worksheet Topics | Mathematical Functions. Кроме "школьных" функций, он включает ряд специальных (например, бесселевых) функций. Такая возможность позволяет делать довольно сложные преобразования сеток.
Бланкирование участков суши. Выше уже неоднократно использовался термин "бланкирование". Объясним, что это такое. На картах водоемов почти всегда есть участки, обозначающие сушу. Это побережье и острова. Эти участки нужно исключить из рассчитанной сетки, поскольку рыбы там, разумеется, нет. Исключение ненужных узлов из расчетов осуществляется присвоением их z-координатам аномального, "дикого" и поэтому легко узнаваемого значения. В программе Surfer это огромное число порядка 1038. Как только в процессе обработки сетки программа натыкается на такое число, она его пропускает, не трогая. Фрагменты суши на карте водоема нередко изображаются чистыми, белыми (blank) участками, отсюда и название "забланкированная область".
Бланкирование осуществляется командой Grid | Blank. После отдачи этой команды нужно указать файл с сеткой, которая будет подвергаться бланкированию, файл границ для участков суши (*.BLN) и файл с сеткой, в который будет записан результат операции бланкирования.
Для того чтобы забланкировать (зачистить) некую группу участков, нужно определить замкнутые контуры их границ. Координаты этих контуров должны храниться в специальном текстовом файле с расширением BLN. Этот файл должен быть подготовлен заранее. Его формат очень прост. Каждый замкнутый участок описывается следующим образом. Сначала идет заголовочная строка,
где даются два числа, а за ними в той же строке — произвольный текстовый комментарий в двойных кавычках. Первое число заголовка определяет количество узловых точек, из которых состоит замкнутая граница участка. Второе число — это код, указывающий программе, куда надо бланкировать: если этот код равен единице, то зачищается все внутри этого участка, если код равен нулю, то зачищается все снаружи участка. В нижеследующих строках идут пары чисел, каждая из которых представляет xy-координаты очередной узловой точки границы. Координаты первой и последней точек должны в точности совпадать, обеспечивая тем самым замыкание границы. В текстовом BLN-файле все замкнутые участки, подлежащие бланкированию, расположены один под другим. Для иллюстрации приведем простой пример такого BLN-файла: 4,1 "low left triangle of shore" 0, 0 0, 3 4, 0 0, 0
5,1 "square island" 4, 5
4, 6
5, 6 5, 5 4, 5.
Здесь описаны две замкнутые области: треугольный участок "материкового побережья" в левом нижнем углу карты и квадратный "остров" в верхней части карты. В результате бланкирования сетки этим BLN-файлом получается карта изолиний, изображенная на рис. 6. Отметим кстати, что незабланкированный вариант этой сетки показан на рис. 1 под другим углом зрения. На рис. 6 особенно ярко выделяются "бычьи глаза" — артефакт метода Inverse Distance to a Power.
Рис. 6. Пример бланкирования сетки треугольным и квадратным участками "суши"
Fig. 6. Example of grid blanking by triangle and square "land" regions
Градусные координаты границ побережий (нулевые изобаты) для всех дальневосточных морей или выбранных участков в этих морях можно получить при помощи специальной программы. Это программа Connect, которая обслуживает базу данных JapOxBer, созданную океанологами ТИНРО-центра.
Бланкировать можно не только участки суши. Допустим, нужно сделать карту распределения биомассы некоего глубоководного вида рыб, о котором известно, что на глубинах менее 500 м он никогда не встречается. Бланкирующей (наружу) границей в таком случае будет изобата 500 м. В принципе нужно бланкировать также и те участки акватории, о которых отсутствует информация, т.е. там, где не было взято ни одной пробы. Интер- и экстраполяция сетки на крупные районы, не обеспеченные фактическим материалом, обычно приводят к непредсказуемым, а поэтому плохим результатам. Выявлять и бланкиро-
вать такие районы можно при помощи некоторых методов расчета сетки (см. обзор выше).
Вычисление объема. Это центральная операция из тех, что нужны для расчета запаса популяции. Несмотря на ее важность, она очень проста. Нужно отдать команду Grid | Volume., указать GRD-файл с построенной сеткой и получить отчет по результатам вычисления объема. Объем пространства, находящегося под сеточной поверхностью, вычисляется сразу тремя методами интегрирования: по правилу трапеций, по правилу Симпсона и по правилу Симпсона 3/8. Обычно все три метода дают совпадающие результаты в первых трех цифрах ответа. Такой точности вполне достаточно. Эта оценка и есть величина промыслового запаса популяции. Здесь же, в отчетном рапорте, в строке Positive Planar Area (Cut) можно узнать общую площадь акватории, по которой оценивался запас.
Масштабирование. Во всех предыдущих рассуждениях предполагалось, что площадь зоны облова для одной пробы и площадь всей акватории, для которой нужно вычислить суммарный запас, выражаются в одних и тех же единицах — например, в квадратных километрах. Если это не так, то оценку суммарного запаса нужно умножить на соответствующий масштабный множитель, приводящий и площадь облова пробы, и площадь акватории к одним и тем же единицам.
Если изучаемая акватория достаточно велика по размерам, например, все море, то теперь уже нужно учитывать шарообразность Земли. Координаты границ у крупных географических объектов выражаются не в километрах, а в градусах. Здесь имеются в виду десятичные градусы, используемые в программе Surfer. Точка, имеющая координату 46 градусов, 30 минут и 00 секунд, в десятичных градусах представлена числом 46,5.
На рис. 7 показана планета Земля с нанесенной на нее градусной сеткой. Эта сетка устроена сложнее, чем простая прямоугольная система координат, предназначенная для небольшого плоского участка земной поверхности. По мере продвижения от экватора к полюсам трапеции, образованные пересечениями параллелей и меридианов, все более сужаются по бокам. Поэтому площадь трапеции с габаритами 1 градус широты на 1 градус долготы, если ее выразить в квадратных километрах, уменьшается от экватора к полюсу. Это обстоятельство надо учитывать при расчете запаса крупных популяций на акватории целых морей. Его нужно также принимать во внимание и для малых акваторий (бухт, заливов, небольших озер), если система координат их береговой линии и координат взятия проб также представлена в градусах, а такая ситуация встречается все чаще в связи с широким распространением спутниковых систем позиционирования.
Рис. 7. Дальневосточные моря России в ортогональной проекции
Fig. 7. Far-Eastern seas of Russia in orthographic projection
Планета Земля — это геоид. Она немного сплющена от полюсов к экватору: ее полярный радиус (расстояние от центра Земли до полюса) на 21 км короче, чем экваториальный (расстояние от центра до экватора). Тем не менее для упрощения расчетов используют идеальную сферу с объемом, равным объему реальной Земли. Эта сфера имеет радиус 6371 км. Неизбежные ошибки, возникающие из-за такой идеализации, при определении расстояний не превышают плюс— минус 0,5 %. Для наших целей такими погрешностями можно пренебречь.
Таким образом, можно принять, что длина дуги размером в один градус широты или долготы равна 111,2 км на экваторе. По мере роста широты и продвижения к полюсам длина такой одноградусной дуги, расположенной на широтном кольце, уменьшается по закону косинуса. При этом длина одноградусной дуги на меридиане остается постоянной, равной тем же самым 111,2 км, т.е. не зависит от удаленности от экватора.
Нам нужно рассчитать среднюю арифметическую между двумя косинусами — для верхней и для нижней широтной параллели у одноградусной трапеции. Вместо этого для упрощения расчетов заменим эту среднюю величину косинусом от средней широты, проходящей через центр этой трапеции. Отличия между этими двумя вариантами весьма малы и начинаются с 4-5-го знака после десятичной точки. Тогда площадь одноградусной трапециевидной ячейки на широтно-меридиональной сетке будет равна 111,2 • 111,2 • cos(a) = 12365 ' cos(a) квадратных километров, где а — широта в градусах, на которой находится центр этой ячейки.
Именно на этот масштабный коэффициент, равный 12365 ' cos(a), надо умножать оценку суммарного запаса популяции. Здесь а теперь уже означает географическую широту, которая проходит примерно по центру изучаемого ареала. Такой масштабный коэффициент необходим в том случае, если площадь зоны облова одной пробы выражается в квадратных километрах, а площадь попу-ляционного ареала — в единицах квадратов градусов. Такие площадные единицы возникают в тех типичных случаях, когда координаты береговой линии моря и местоположения пробных уловов (используемые программой Surfer для расчета сетки) представлены в десятичных градусах. Дополнительное изменение величины предлагаемого здесь масштабного коэффициента при площадях пробы, выражаемых не в квадратных километрах, а в гектарах, квадратных метрах и т.д., обсуждать излишне.
Точно такой же масштабный коэффициент, равный 12365 ' cos(a), следует использовать для пересчета площади ареала из квадратов градусов в квадратные километры. Об оценивании площади ареала говорилось выше. Оценка находится в строке Positive Planar Area (Cut) в отчетном рапорте по расчету площадей и объемов.
Пример
Важное, но редко упоминаемое ограничение для применимости площадного метода Аксютиной и его всех средневзвешенных модификаций состоит в следующем. И размер пробы, и размер всей акватории должны выражаться в площадных единицах. Это означает, что при расчетах можно вообще не учитывать закономерности в вертикальном распределении особей по глубине. Если же такой учет и придется делать, то это должно проходить просто и легко. Такое условие применимо для популяций бентосных или придонных видов, а также для популяций, обитающих в достаточно тонком слое воды, который по вертикали полностью или в значительной степени охватывается орудием лова. В данной статье приводится пример расчета запаса для популяции рыб, подчиняющейся этому ограничению.
Материал. Выбранный для иллюстрации вид периодически заходит в российскую экономическую зону Японского моря. Это дальневосточная сардина Sardinops melanostictus.
Используем материал, собранный сотрудниками лаборатории прикладной биоценологии ТИНРО-центра. С 1981 по 2003 г. в российской и частично корейской экономической зоне Японского моря было проведено 32 научно-исследовательских рейса. Не весь собранный материал годится для оценивания запаса именно дальневосточной сардины. Во-первых, из исходных данных нужно исключить сведения по последним годам наблюдений: сардина ловилась в российской экономической зоне Японского моря лишь по 1990 г. включительно. Во-вторых, не нужно использовать данные, полученные с декабря по апрель: нагуливающаяся сардина обитает в районе исследований лишь в теплое время года. В-третьих, следует исключить траления, проводившиеся в темное время суток (с 18 часов вечера до 6 утра): сардина ночью практически не ловится. И наконец, в-четвертых, не надо использовать те пробы, средняя глубина траления которых превышала 100-метровую отметку, чтобы исключить из анализа косые траления, в которые сардина попадала лишь на последних этапах подъема трала. После такой многоступенчатой селекции материала из исходных 2119 стандартных часовых тралений в нашем распоряжении остается лишь 667. При определении биомасс был учтен коэффициент уловистости, для сардины принятый равным 0,4 (Нектон ..., 2004).
Вероятностное распределение япономорской сардины-иваси по массе тела подчиняется логнормальному закону с "сигмой" 0,869 и с модой 148 г. Проинтегрировав плотность логнормального распределения с этими параметрами, нетрудно найти долю особей, масса которых превышает заданный промысловый порог. Таким образом, если оценку запаса всей сардины в исследованном регионе умножить на эту долю, то легко получить оценку промыслового запаса. Тем не менее такую оценку мы не проводили. Промысловый размер сардины для рыбаков не устанавливался вообще, поскольку ими облавливалась только нагульная часть популяции с крупнотелыми особями.
Представленный пример является химерой (".перед львиный, средина козлиная, а зад змеиный."), поскольку здесь в одном массиве объединена разномасштабная временная изменчивость в пространственном распределении обилия этого вида. Но для иллюстративных целей методической статьи этот пример годится.
Распределение по глубине. Это важная экологическая характеристика вида, оказывающая существенное влияние на оценку запаса. Нередко вертикальный размер орудия лова не охватывает весь диапазон глубин, в котором обитает данный вид. В таких случаях требуется расчет специальных поправок на улови-стость по вертикали. К обсуждению этих поправок мы сейчас и перейдем.
Нужно признаться, что именно эта часть работы оказалась самой трудной. В процессе ее выполнения автору удалось совершить практически все серьезные ошибки, какие здесь вообще можно придумать. В первом варианте было построено обычное вероятностное распределение для биомассы сардины по глубине: на оси абсцисс — средняя глубина траления, на оси ординат — доля биомассы в том или ином диапазоне глубин. Потом стало понятно, что всю пойманную данным конкретным тралом биомассу нельзя сосредотачивать исключительно на той глубине, где находится центр устья трала. Нужно хотя бы равномерно распределить данную биомассу по всему вертикальному раскрытию этого трала. И наконец, пришлось учесть то обстоятельство, что частота распределения самих тралений практически всегда неодинакова по глубине — безотносительно от того, сколько биомассы попадает в тралы на разных глубинах (рис. 8).
Типичное унимодальное частотное распределение самих тралений порождало бы абсолютно такое же распределение биомассы по глубине, если бы в каждый трал попадало в точности одно и то же постоянное количество рыб. В реальности количество пойманных рыб и их биомасса варьируют от трала к тралу. Поэтому вероятностное распределение биомассы по глубине имеет другие пока-
затели центра и разброса. Тем не менее первое, что приходит в голову — местоположение моды (наибольшей доли) на оси глубин в этом распределении биомасс указывает на некую оптимальную глубину обитания вида. На самом деле этот "оптимум" обычно имеет весьма отдаленное отношение к экологической реальности. Проиллюстрируем сказанное образчиком классического английского юмора.
300
Рис. 8. Частотное распределение числа тралений по глубине верхней подборы трала. Учтены только те тралы, которые удовлетворяют принятым для сардины условиям отбора материала
Fig. 8. Frequency distribution of creeping number along the depth of upper drag-net boundary. Creepings are taken into account only when conditions of data selection for sardine are satisfied
>s в в
OJ
& H
о ч о В =T
200
100
5 10 15 20 25 30 35 40 Глубина, м
Геоботаник П. Грейг-Смит (1967) приводит пример, позаимствованный у Эшби (Ashby, 1936), который показал, как найти оптимальную величину рН почвы для произрастания телеграфных столбов. Для этого следует взять побольше проб почвы под столбами и построить частотное распределение этих проб по величине рН. Значение рН, приходящееся на максимальную частоту встречаемости, указывает на искомый оптимум кислотности почвы, наиболее благоприятный для произрастания телеграфных столбов.
Для того чтобы избавиться от искажающего влияния частоты тралений на разных глубинах, при построении вертикального распределения биомассы нужно организовать так называемую равномерную выборку, в которой каждый класс глубин должен быть представлен одинаковым числом тралений. Выбор конкретной пробы из полного списка тралений, попавших в данный диапазон глубин, можно делать либо регулярным способом, последовательно продвигаясь вниз по этому списку (может быть, и неоднократно, чтобы соблюсти равноправность разных классов глубин), либо организовав из этого списка случайный выбор с возвращением при помощи генератора случайных чисел. Здесь использован последний способ.
На рис. 9 приведено вероятностное распределение биомассы уловов сардины по глубине в северо-западной части Японского моря. Анализ этого распределения показывает, что оно ограничено слоем верхней эпипелагиали от 0 до 7080 м. Средняя глубина обитания сардины равна 32,1 м, а среднеквадратическое отклонение у этого распределения равно 14,1 м.
Рис. 9. Нормированное распределение биомассы сардины по глубине. Сумма всех долей равна единице. Плавная кривая — модель Вейбулла (1)
Fig. 9. Normalized depth distribution of sardine biomass. Sum of all fractions is equal to unity. Smoothed line is the Weibull model (1)
20 40
Глубина, м
0
0
Для аппроксимации распределения биомассы по глубине x была выбрана модель Вейбулла (Хастингс, Пикок, 1980):
f (x) = (cxc-1/bc) • exp [~(x/b)) ] • dx, (1)
где b, c — параметры модели, dx — классовый промежуток плотности распределения, в нашем примере равный 5 м.
Была сделана подгонка модели Вейбулла f(x) под плотность распределения биомасс по глубине x. Эта процедура легко осуществляется во многих статистических программах. Оценки модельных параметров: b = 35,58 ± 0,44 м, c = 2,607 ± ± 0,066, коэффициент детерминации модели достаточно велик и равен 0,98.
В этой модели можно учитывать суточные вертикальные миграции популяции. Для этого надо оценить модельные параметры b и c отдельно для разных часовых интервалов, а затем аппроксимировать суточную динамику коэффициентов b и c при помощи синусоид. Время суток изменяется, синусоиды для коэффициентов b и c в модели (1) колеблются, поэтому распределение биомассы по глубине демонстрирует циклические колебания с суточным периодом. Точно такой же подход можно использовать для описания сезонных вертикальных миграций. В нашем конкретном примере такое развитие модели не понадобилось: у нагуливающейся дальневосточной сардины нет заметных вертикальных миграций.
Кумулята (огива) распределения Вейбулла легко вычисляется как
F (x) = 1 - exp [-(x/b)) ]. (2)
Величина F(x) есть доля от всей биомассы, содержащаяся в диапазоне глубин от 0 до x метров.
Доля биомассы, захваченной тралом из всего ее распределения по глубине, вычисляется двумя подстановками в кумуляту (2). Сначала подставляем туда вместо x глубину нижней подборы трала xmax. Получаем долю Fmax. Затем в ту же формулу кумуляты вместо x подставляем глубину верхней подборы xmin. Получаем долю F . . Тогда доля захваченной тралом биомассы равна P = F - F . . На
min 1 1 max min
рис. 10 показан пример с аргументами xmin = 20 м, xmax = 50 м. При найденных ранее параметрах b = 35,58 ± 0,44 и c = 2,607 ± 0,066 получаем оценку
P = exp [-(xmin/ b) ]-exp [-(xmffi/b )c ] = 0,800 - 0,088 = 0,712. Иными словами, в данном примере расчета трал с вертикальным размахом устья от 20 до 50 м глубины захватывает чуть более 2/3 от всей биомассы сардины в этом месте.
Рис. 10. Определение доли Р захваченной тралом биомассы сардины из всего ее распределения по глубине: xmin — глубина верхней подборы трала, x — глубина нижней
1 1 max J
подборы
Fig. 10. Determination of catched fraction P of sardine biomass from all its depth distribution: x . is upper
1 min 1 1
boundary of trawl mouth, x is lower
max
boundary of trawl mouth Глубина, м
Статистическая ошибка этой оценки может быть вычислена по простой, хотя и громоздкой формуле переноса ошибок (Худсон, 1970) — особенно простой, если учесть, что корреляция между ошибками модельных параметров b и c
оказалась пренебрежимо мала (0,009). Мы не будем здесь отвлекаться на эту второстепенную тему и сразу дадим результат: P = 0,712 ± 0,012.
Теперь у нас есть все сведения, необходимые для расчета. Пусть, например, в данном улове оказалось поймано 10 т сардины. Тогда ее полная биомасса в этом месте должна быть равной (10/0,712) ± (10 • 0,012/0,7122) = 14,04 + 0,24 т по всему диапазону глубин. Менее чем 2 %-ная ошибка этой оценки типична для обсуждаемого материала по сардине, и далее ее можно совсем не учитывать в расчетах.
По описанной методике были рассчитаны полные биомассы для каждого из
667 тралений с его уникальными значениями биомассы улова и глубинами x . , r J / \ т xmax. При этом величина xmax трактовалась как эффективная (приведенная) глубина нижней подборы. Это была нижняя горизонтальная граница прямоугольника, ширина которого принималась равной горизонтальному раскрытию трала, а площадь приравнивалась к площади реального эллипсовидного устья трала. В среднем по всем тралениям доля вертикального захвата биомассы тралом составила для сардины P = 0,638.
Расчет поверхности уловов. Рассчитанные по описанной выше методике данные по биомассам сардины пошли на создание карты распределения этого показателя по акватории района. Выбор алгоритма расчета сетки сначала осуществлялся на основе метода перекрестной проверки (Cross Validation — см. выше). Этот метод выделил группу лидеров с наименьшими значениями ошибки аппроксимации. Все алгоритмы, попавшие в группу лидеров, мало отличались друг от друга по этому критерию отбора. Поэтому окончательный выбор среди них наилучшего алгоритма был сделан из других соображений, основанных на рекомендациях разработчиков программы Surfer.
Был выбран метод кригинг с шагом сетки 0,1о и с последующей 50-кратной ее фильтрацией при помощи стандартного гауссовского фильтра. Результат показан на рис. 11. Это ареал дальневосточной сардины на российской и частично корейской акватории, рассчи
танный на основе полных биомасс уловов по всему диапазону глубин в верхней эпипе-лагиали.
Рис. 11. Распределение биомассы сардины (т/км2) в северо-западной части Японского моря. Точки — места проведения тралений
Fig. 11. Sardine biomass distribution (ton/km2) in the northwestern of Japan Sea. Points are the locations of trawlings
52 51504948 47 46 4544 43 42 41 40 39
. Sardinops. m melanostiçtus
1 300
— 150
— 80
— 40
— 20
— 10
— 5
— 2
0
■ te * I
■.......f
: '
.........ii
128 130
151
132 134 136 138 140 142
Запас. Центр исследованного ареала сардины проходит по широте около 45,25°. Отсюда масштабный коэффициент равен 12365 ' cos(45,25°) = 8705 км2 град-2. Площадь изученного ареала составляет 44,77 град2. Умножая эту величину на 8705 км2 град-2, получаем, что площадь обследованной акватории равна 390 тыс. км2.
Общий объем, накрываемый сеточной поверхностью биомасс, оказался равным 2562 т ' км-2 ' град2. Умножим эту оценку на масштабный коэффициент 8705 км2 град-2 и получим запас сардины, равный 22,30 млн т. Среднеквадрати-ческую ошибку средней, а значит и ошибку оценки запаса, здесь вычислить невозможно, так как количество узлов в рассчитываемой сетке задается пользователем, и его нельзя трактовать как объем выборки N.
Средний улов по "сырым", первичным исходным данным составляет 47,19 ± ± 9,19 т ' км-2. Умножая эту оценку на площадь акватории (390 тыс. км2), получаем запас, вычисленный по классическому площадному методу Аксютиной. Он равен 18,40 ± 3,58 млн т.
Степенное преобразование данных. Нужно отметить, что частотное распределение уловов сардины имеет L-образный вид и характеризуется резко выдающимся правым "хвостом": чем крупнее улов, тем меньше шансов его получить. Известно, что чем длиннее этот хвост, тем сильнее повышается оценка запаса популяции. Для ослабления этого эффекта З.М. Аксютина (1968) предлагает вычислять не среднюю арифметическую по уловам, а среднюю геометрическую. Последний параметр оценивается следующим образом: средняя считается не по исходным уловам, а по их логарифмам, после чего этот средний логарифм потенцируется. Такой прием действительно ослабляет завышающий эффект правого хвоста. К сожалению, предложенное З.М. Аксютиной логарифмическое преобразование обладает серьезным недостатком: логарифм нельзя брать от нуля. Вместе с тем нулевые уловы несут в себе важную информацию об ареале популяции (в бразильской сельве не встретишь чукотский чум — очень жарко, однако). Поэтому заменять нули произвольно заданной малой константой или отбрасывать их вообще — это, очевидно, неправильно.
Вместо логарифмического можно использовать степенное преобразование вида u = zk со степенью меньше единицы (k < 1). Оно также выпукло вверх и поэтому ослабляет эффект правого хвоста в частотном распределении биомасс. Кроме того, примененное к нулю такое преобразование оставляет этот нуль без изменений, а не порождает минус-бесконечность.
По рекомендации С.А. Пионтковского (2003) и Даунинга с соавторами (Downing et al., 1987), конкретное значение степени в степенном преобразовании уловов следует взять равным 0,2. Тогда дисперсия перестает зависеть от условной средней, и материал приобретает полезное свойство гомоскедастично-сти. Мы сочли достаточной степень, равную 1/4. Она почти такая же по величине, как 0,2, но в электронных таблицах легко реализуется как двойной квадратный корень sqrt(sqrt(...)), не требуя логарифмирования и последующего потенцирования.
Был повторен расчет поверхности, описывающей пространственное распределение биомасс сардины. Но на этот раз использовались нелинейные преобразования исходного материала. Вначале была получена сетка, описывающая пространственное распределение корней четвертой степени из биомасс сардины. Затем, после фильтрации и бланкирования, эта сетка была подвергнута обратному степенному преобразованию: каждое z-значение узла возводилось в четвертую степень, т.е. четырехкратно умножалось само на себя.
Суммарный объем, накрываемый сеточной поверхностью нелинейно преобразованных биомасс, на этот раз оказался равным 597,3 т ' км-2 ' град2. Умножая эту оценку на масштабный коэффициент 8705 км2 град-2, получаем, что запас сардины, рассчитанный по нашей методике, составляет 5,20 млн т. Эта оценка
152
вчетверо ниже оценки запаса в 22,30 млн т, которая была получена ранее по биомассам, не подвергнутым прямому и обратному степенному преобразованию.
Первичные исходные материалы по уловам также можно подвергнуть степенному преобразованию со степенью 1/4, на основе этих данных вычислить среднюю, которую затем подвергнуть обратному степенному преобразованию, возведя в четвертую степень. Рассчитанный таким способом средний улов оказывается равным 2,025 т ' км-2 . Умножая эту величину на площадь акватории, получаем запас сардины по модифицированному площадному методу Аксюти-ной. Он равен всего лишь 0,79 млн т.
Сравнение с независимой оценкой. Можно сопоставить наши расчеты с независимой информацией. В работе В.А. Дударева (1985) приводится величина запаса сардины в 1982 г., основанная на применении акустических методов. Этот запас равен 6,7 млн т. Здесь нужно учесть, что биомасса сардины, усредненная за весь использованный нами период наблюдений с 1981 по 1990 г., составляла лишь 72,3 % от биомассы за 1982 г., поэтому оценку В.А. Дударева надо скорректировать с учетом этого факта. Получаем усредненный запас, равный 6,7 ' 0,723 = 4,85 млн т.
Оценка, подсчитанная по классическому площадному методу З.М. Аксюти-ной (18,4 млн т), в 3,8 раза превысила запас по В.А. Дудареву. Оценка по модифицированному методу З.М. Аксютиной (0,79 млн т), напротив, оказалась в 6,1 раза ниже запаса по В.А. Дудареву.
Запас, подсчитанный по объему под поверхностью уловов, не подвергнутых степенному преобразованию (22,3 млн т), оказался в 4,6 раза выше, чем запас по В.А. Дудареву.
И наконец, ближе всего к запасу по В.А. Дудареву (4,85 млн т), найденному акустическими методами, оказалась наша оценка запаса сардины (5,20 млн т), полученная расчетом объема под поверхностью уловов, которые сначала были трансформированы прямым, а затем после фильтрации — обратным степенным преобразованием. Различия между этими двумя оценками составляют всего лишь 7 %.
Мы никогда не узнаем истинной величины запаса сардины в северо-западной части Японского моря. Тем не менее для того чтобы отбирать разные методы, дающие близкие результаты, подобные сравнения и взаимные проверки нужно делать постоянно на разных видах и в разных водоемах. Пока же из четырех опробованных здесь способов расчета запаса наилучший результат показал метод расчета объема под поверхностью, которая описывает пространственное распределение уловов, подвергнутых степенному преобразованию.
Автор выражает признательность главному научному сотруднику ТИНРО-центра д-ру биол. наук, профессору В.П. Шунтову — за предложение написать настоящую статью (сам я как-то на это не решался); ведущему научному сотруднику лаборатории прикладной биоценологии ТИНРО-центра, канд. биол. наук О.А. Иванову — за возможность проделать нужные расчеты при помощи лицензионной копии программы Surfer ver. 8; всем сотрудникам лаборатории прикладной биоценологии ТИНРО-центра — за материал по дальневосточной сардине Японского моря.
Список литературы
Аксютина З.М. Элементы математической оценки результатов наблюдений в биологических и рыбохозяйственных исследованиях. — М.: Пищ. пром-сть, 1968. — 289 с.
Борисовец Е'Э', Вдовин А.Н., Панченко В.В. Оценки запасов керчаков по данным учетных траловых съемок залива Петра Великого // Вопр. рыб-ва. — 2003. — Т. 4, № 1. — С. 157-170.
Грейг-Смит П. Количественная экология растений. — М.: Мир, 1967. — 359 с.
Дударев В.А. Сардина Японского моря, ее экология и промысел: Дис... канд. биол. наук. — Владивосток: ТИНРО, 1985. — 179 с.
Мальцев В.А. К вопросу о получении несмещенных оценок при моделировании месторождений с нелинейными преобразованиями данных // Горный информационно-аналитический бюллетень: Ершовские чтения. — 1999. — Вып. 5. — С. 106-111.
Нектон северо-западной части Японского моря. Таблицы численности, биомассы и соотношения видов / под ред. В.П. Шунтова и Л.Н. Бочарова. — Владивосток: ТИНРО-центр, 2004. — 226 с.
Никольский Г.В. Теория динамики стада рыб, как биологическая основа рациональной эксплуатации и воспроизводства рыбных ресурсов. — М.: Наука, 1965. — 383 с.
Пионтковский С.А. Многомасштабная изменчивость мезопланктонных полей океана: Дис.д-ра биол. наук. — Севастополь, 2003. — 326 с.
Столяренко Д.А., Иванов Б.Г. Метод сплайн-аппроксимации плотности для оценки запасов по результатам траловых донных съемок на примере креветки Pandalus borealis у Шпицбергена // Морские промысловые беспозвоночные. — М.: ВНИРО, 1988. — С. 45-70.
Суханов В.В. Научная графика на компьютере. — Владивосток: Дальнаука, 2005. — 355 с.
Хастингс Н., Пикок Дж. Справочник по статистическим распределениям. — М.: Статистика, 1980. — 95 с.
Худсон Д. Статистика для физиков. — М.: Мир, 1970. — 296 с.
Ashby E. Statistical ecology // Bot. Rev. — 1936. — Vol. 2. — P. 221-235.
Downing J.A., Perusse M., Preuette J. Effect of interreplicate variance on zooplankton sampling design and data analysis // Limnol. Oceanogr. — 1987. — Vol. 32, № 3. — P. 673-680.
Quenouille M.H. Notes on bias in estimation // Biometrica. — 1956. — Vol. 43. — P. 353-360.
Ripley B.D. Spacial statistics. — N.Y.: Wiley, 1981.
Sibson R. A vector identity for the Dirichlet Tesselation // Math. Proc. Cambridge Phil. Soc. — 1980. — Vol. 87. — P. 151-155.
Поступила в редакцию 17.12.07 г.
