DOI: https://doi.org/10.37162/2618-9631-2022-3-42-77

УДК 551.509.313+551.509.324.2+551.508.85

Верификация радиолокационного наукастинга
областей осадков значительной площади
с помощью обобщенного распределения Парето.
Часть 2: приложение к прогнозам в теплый
и холодный периоды 2017–2018 гг.

А.В. Муравьев, А.Ю. Бундель, Д.Б. Киктев, А.В. Смирнов

Гидрометеорологический научно-исследовательский центр
Российской Федерации, г. Москва, Россия

muravev@mecom.ru

Обобщенное распределение Парето использовано для моделирования распределения значений площадей тех осадков, которые наблюдались и прогнозировались по зонам обзора отдельных радиолокаторов схемой наукастинга Гидрометцентра России в теплый и холодный периоды 2017–2018 гг. Протестированы различные методы оценок параметров распределения и доверительных интервалов. Основное внимание уделено параметру формы, определяющему поведение хвостов распределений. Обобщенная оценка качества прогноза построена на доли пересечения соответствующих доверительных интервалов. Показано, что для большинства случаев порог Парето в 625 точек, равносильный квадрату 50×50 км, отделяет такие объекты большего размера, которые удовлетворительно моделируются тяжелохвостым распределением и которые по обоим периодам года вполне качественно (в "климатологическом" среднем) прогнозируются системой наукастинга осадков.

Ключевые слова: наукастинг осадков, пространственная верификация, оценки осадков радиолокационными методами, статистический анализ превышения порогов, математические пакеты для анализа экстремальных величин

Verification of radar precipitation nowcasting
of significant areas
using generalized Pareto distribution.
Part 2: application to forecasts in warm
and cold periods of 2017–2018

A.V. Muravev, A.Yu. Bundel, D.B. Kiktev, A.V. Smirnov

Hydrometeorological Research Center of Russian Federation, Moscow, Russia

muravev@mecom.ru

The generalized Pareto distribution was used to model the distribution of precipitation area sizes that were observed and predicted for the coverage zones of individual radars by the Hydrometeorological Research Center of the Russian Federation's nowcasting scheme in 2017-2018. Various methods for estimating the distribution parameters and confidence intervals were tested. The main attention was paid to the estimates of the shape parameter that determines the behavior of the distribution tail. The generalized assessment of the nowcasting quality is built on the intersection ratio of the corresponding confidence intervals. It is shown that for most cases the Pareto threshold of 625 points, which is equivalent to a square of 50´50 km, separates objects of larger sizes that are satisfactorily modeled by the heavy-tailed distribution and which are quite acceptably (on "climatological" average) predicted by the precipitation nowcasting system for specific periods of the year.

Keywords: precipitation nowcasting, spatial verification, radar precipitation estimates, statistical analysis of threshold exceedances, mathematical packages for extreme value analysis

Введение

В представленной статье модель "пиков над порогом", описываемая обобщенным распределением Парето и подробно обсужденная в первой части статьи [5], применяется к значительным и максимальным размерам областей осадков.

Исходными данными послужили архивы испытаний системы радиолокационного наукастинга осадков в теплый (май–сентябрь 2017 г.) и в холодный (ноябрь 2017 – март 2018 г.) периоды года. Данные сгруппированы по обзорам отдельных радиолокаторов ДМРЛ-С на территории Центрального федерального округа. По каждому периоду накоплено порядка 20 тысяч полей наблюдений и прогнозов на 120 минут в 10-минутной дискретности, при этом все поля, содержащие вписанные круги радиолокационных обзоров с радиусом около 250 км, реорганизованы в матрицы размером 256´256.

Для приложения методов объектно-ориентированной верификации области осадков (вообще говоря, разрывные и с фрактальными границами) преобразуются в односвязные объекты (области без внутренних "дырок") с помощью подходящей изогиеты, пространственного осреднения и некоторых приемов слияния изолированных областей малого масштаба. Топологические понятия связности и односвязности областей евклидова пространства вполне покрываются для наших целей словосочетанием "область сплошных осадков". Определение, идентификация и характеризация объектов производятся с помощью пакета SpatialVx в среде языка R.

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

Модель пиков над порогом более полувека используется при анализе экстремумов в гидрологических задачах, в которых оцениваются экстремальные характеристики слоев осадков, их площадей и длительности ливней. Ссылки на тематические руководства ВМО по гидрологии также приведены в [5].

Некоторых пояснений требует оценка площади области через подсчет точек внутри и на границе области. При таком подходе случайной величиной оказывается неотрицательное целое число, и распределение становится по существу дискретным. Использование такого "размера" в качестве оценки площади допустимо благодаря эквидистантной сетке радиолокационного поля и относительно небольшой погрешности, вносимой граничными точками выделяемой области. Имеется достаточное количество алгоритмов для расчета фактической площади областей на земной поверхности (например, [4]), но мы довольствуемся этим простейшим способом, который, кстати, плодотворно применялся в гидрологических исследованиях до эпохи интернета (например, [15, 16]). Ниже термины размер и площадь объекта будут использоваться как эквивалентные в указанном выше смысле.

Само дискретное распределение Парето обнаруживается, часто под названием закона Ципфа, в самых разных явлениях, помимо рассмотренного В. Парето распределения доходов населения выше некоторого уровня. Сюда также относятся распределения городов по количеству жителей, биологических родов по количеству видов, ученых по количеству публикаций, слов в произведениях по частоте употребления. К середине прошлого века была построена модель "случайного роста" [14], порождающая все перечисленные в этом абзаце дискретные распределения единым алгоритмом, сходным с методом построения марковской цепи с помощью приемов Монте-Карло (МСМС). Построенный временной ряд случайного роста является вариантом случайного блуждания под названием модель Саймона –Юла.

Дополнительно прокомментируем в терминах теории верификации использованную нами оценку прогнозов осадков с помощью распределения существенных (максимальных) площадей областей осадков. Речь идет об оценке качества прогнозов с привлечением не полной двумерной плотности распределения (по методологии Мерфи – Винклера), а лишь "маргинальных" распределений, или "выборочных климатологий", рассчитанных отдельно для наблюдений и прогнозов. При этом вначале оцениваются параметры распределения Парето, а затем попарно сравниваются между собой оценки полученных параметров масштаба и формы. Из известных атрибутов верификации прогнозов восстановимы "точность" (accuracy, сравнение средних по обеим выборкам прогнозов и наблюдений), "смещение" (bias, разность средних), "категоричность" (sharpness, дисперсия по прогнозам) и "неопределенность" (uncertainty, дисперсия по наблюдениям). Средние и дисперсии однозначно восстанавливаются по оценке параметров Парето с единственным исключением случаев тяжелого хвоста, когда нужных моментов распределения не существует; выборочные оценки моментов можно рассчитать, но они будут некорректны. Совершенно очевидно, что наилучшее качество наукастинга областей осадков значительной площади достигается в тех случаях, когда такого рода области уже сформированы и включены в начальные поля статистической схемы наукастинга: схема генерирует небольшие отклонения от среднего лишь в блоке стохастического возмущения, имитирующего стандартные допуски для оценок радиолокационной отражаемости. Опора на инерцию метеорологического процесса вполне оправдана, так как, во-первых, на сроках до двух-трех часов гидродинамические модели еще не могут предложить продукцию конкурентного качества, во-вторых, статистическая схема радиолокационного наукастинга проста, экономична и обновляет, а значит, и уточняет свой собственный прогноз каждые десять минут.

Нет нужды специально подчеркивать, что современное статистическое моделирование вообще и экстремумов в частности считается обоснованным и достоверным только при опоре на математические пакеты, прошедшие экспертизу профессионального сообщества. Использованные нами расчетные пакеты SpatialVx и extRemes были созданы для среды языка R, обстоятельно документированы, а результаты приложений опубликованы в солидных журналах. Их важнейшие характеристики обсуждены нами в [5], здесь будут добавлены сведения по заданию некоторых опций пакета extRemes [11].

1. Настройка опций для функций расчетного пакета extRemes

Общие характеристики. Детальное описание extRemes и учебные материалы можно найти в [9–11]. Данный пакет (на жаргоне языка R – библиотека) содержит разнообразные функции для анализа экстремумов с возможностью подключения дополнительных переменных (ковариат) и с приемами декластеризации интервалов по методологии [6].

В нашей работе использованы две функции пакета, fevd() и ci(), и две вспомогательные функции языка R для выдачи выходной продукции, summary() и print(). Функциональные скобки указывают на наличие настроечных опций.

Перечислим названия четырех методов оценивания параметров в функции fevd, детально описанные в первой части [5]. Дадим их расшифровку, часто оставляя сокращения в тексте для экономии места: метод максимального правдоподобия (MLE), метод L-моментов (L-moments), "обобщенный" метод максимального правдоподобия (GMLE) и байесовский метод (Bayesian). Статистический вывод (расчет доверительных и достоверных интервалов) проводится с использованием 1) нормального распределения, 2) профиля правдоподобия, 3) байесовских оценок по апостериорной функции, 4) бутстрепа. Функция fevd() имеет опции двумерного анализа экстремальных зависимостей по рекомендациям [6, 13] и дополнительной проверки экстремальности по методике [7].

Оценка качества модели (подгонки параметров) рассчитывается с помощью критерия Акаике (AIC), Байесовского информационного критерия (BIC) и информационного критерия отклонения (DIC) (детали вычислений критериев в [8]). Однако функция fevd не рассчитывает единый критерий для всех методов оценки параметров, ввиду чего нами использован общий критерий Хи-квадрат в приложении к гистограммам, которые строятся базовой функцией hist() языка R.

Задание начальных значений и основных опций функции fevd(). В случае методов MLE/GMLE желательно иметь хорошие начальные оценки параметров, но по умолчанию функция fevd пытается сама отыскать эти приближения по стандартным общим рекомендациям (например, по [13]), перечисленным в первой части нашей статьи.

По умолчанию вначале рассчитываются оценки L-моментов и оценки, основанные на моментах распределения Гумбеля; те оценки, на которых логарифмическое правдоподобие принимает наименьшее (отрицательное) значение, используются дальше. Рекомендуется воспользоваться методом Bayesian и протестировать несколько начальных значений, чтобы удостовериться в отсутствии их влияния на конечный результат; если тестирование не удовлетворяет, прибегнуть к проверенному средству – методу максимального правдоподобия (MLE). Оценка интенсивности потока Пуассона для обобщенного метода правдоподобия (GMLE) рассчитывается по размеру полной выборки (заданием общего объема в 20000).

Так как в Монте-Карловской имитации цепи Маркова имеется вредный "спинап", первые несколько сотен байесовских реализаций следует удалить из результирующей выборки с помощью опции отжига дефекта (burn-in). По умолчанию количество шагов в процессе случайного блуждания равно 10000 с отбраковкой 500 первых реализаций.

Для учета временных и календарных характеристик имеются опции, задающие базовый период (по умолчанию “years”) и единицы времени (по умолчанию “days”).

Для методов GMLE и Bayesian есть возможность заказать априорную функцию распределения параметра формы. По умолчанию для GMLE используется бета-распределение на интервале от -0.5 до 0.5 с параметрами самого распределения p=9 и q=6 (при которых математическое ожидание смещено в положительную область и равно 0.1), а в байесовском методе по умолчанию используются нормальные функции распределения с оценками средних методом MLЕ со стандартным отклонением 10 для всех параметров. Только для оценки методом Bayesian имеется символьная переменная, заказывающая название функции, генерирующей параметры предложения (предложения кандидата) на каждой итерации MCMC. По умолчанию используется цепь случайного блуждания, и при текущем значении параметра предлагается кандидат в виде аддитивной добавки к текущему значению числа, извлеченной из нормального распределения с выборочными оценками ожидания и дисперсии.

Оценка доверительных интервалов. Функция ci() рассчитывает по выходной продукции fevd() доверительные и достоверные (в байесовской оценке) интервалы (ДИ), учитывая метод оценки параметров в fevd() и указанный в ci() метод расчета интервалов.

Для L-moments единственный доступный в пакете метод – параметрический бутстреп с количеством итераций (реплик) R. Рекомендуется определять значение опции R через пробы и ошибки, скажем, начав с R=100, и постепенно увеличивать ее (на одну-две сотни), пока результаты не стабилизируются. По умолчанию R=100.

Для MLE/GMLE можно задать аппроксимацию нормальным распределением. Для расчета ДИ оценок параметров методом Bayesian используются крайние процентили из результирующей выборки MCMC после удаления первых значений.

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

Выходная цифровая и графическая продукция. Оценки параметров и их стандартные ошибки размещаются в векторе выходной продукции функции fevd. В отдельном массиве содержится матрица кросс-ковариаций параметров, а оценки качества аппроксимации доступны в переменных с аналогичными наименованиями (AIC, BIC и DIC).

Характеристики выходной информации задаются опциями базовых функций summary и print языка R. При байесовской оценке создается панель графиков с указанием апостериорной функции распределения и панель следов метода MCMC для каждого параметра. По умолчанию на апостериорной выборке в качестве оценки параметра рассчитывается результирующее среднее, но возможно задание опции для выведения на график апостериорной моды выборки.

2. Оценки качества наукастинга с помощью параметров распределения Парето

Функция fevd() вместе с оценками параметров масштаба и формы возвращает оценки стандартной ошибки, которые ниже будут использованы для расчета границ доверительных интервалов (L,U) – основного инструмента обобщенной оценки качества. Вначале опишем выходную табличную и графическую продукцию, накопленную в результате испытаний системы наукастинга осадков в теплый и холодный периоды 2017–2018 гг.

2.1. Иллюстративные результаты верификации наукастинга областей осадков

Для каждого из восьми радиолокаторов, каждого из двух периодов года, каждого из четырех порогов Парето формируются по две таблицы числовых данных (размеры максимальных объектов и гистограммы) и по одному набору графиков.

Выходные табличные данные. Каждая таблица с характеристиками площадей состоит из шести столбцов: 1) названия, 2) наблюдения на локаторе, 3) – 6) прогнозы по модели STEPS на 30, 60, 90 и 120 мин.

По строкам выделяются две группы характеристик: 1) общие характеристики площадей объектов и для четырех методов оценок параметров, 2) критерии качества моделирования (Акаике и Хи-квадрат) и оценки параметров масштаба и формы с границами их доверительных интервалов.

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

Общие характеристики площадей максимальных осадков. Перейдем к иллюстративным и характерным примерам. Табл. 1 представляет общие статистические свойства площадей объектов по данным локатора RUDL (Смоленск), размером не менее 625 точек, для прогнозов на четыре срока.

Таблица 1. Общие характеристики максимальных площадей объектов, содержащихся в непрерывных ситуациях полей осадков по наблюдениям (obs) ДМРЛ-С Смоленск (RUDL) в теплый период года (май-сентябрь 2017 г.) и в прогностических полях осадков (30-120 мин) по модели STEPS

Table 1. General characteristics of the areas of maximum objects in continuous precipitation situations according to observations (obs) of DMRL-S radar of Smolensk (RUDL) in the warm season (May-September 2017) and in STEPS nowcasting precipitation fields (30-120 min)

	obs 30 60 90 120
peaks	92 119 123 131 130
area_min	625 654 630 630 639
area_med	2496 2394 2380 2195 2150
area_max	19812 19982 20023 19579 18426
ndegf	11 11 12 11 11

Примечание. Названия строк: peaks – количество максимальных объектов (или ситуаций); area_min/ med / max – минимальная, медианная и максимальная площадь, соответственно; ndegf – количество градаций на соответствующих гистограммах (или количество степеней свободы для критерия Хи-квадрат "по правилу Стурджеса").

Характеристики табл. 1 отражают общие свойства выделяемых объектов значительной площади, присущие зонам обзора всех радиолокаторов (за редким исключением).

Количество максимальных по площади объектов в прогностических полях превышает аналогичное количество в полях наблюдений, примерно на 20–25 %. Но в таких статистиках, как среднее, медианное и максимальное значение этих площадей, наборы полей наблюдений и полей прогнозов вполне сопоставимы.

В данных зоны обзора локатора Смоленск наблюдаются следующие особенности. Медианная площадь в полях наблюдений составляет 2496 точек, а в прогностических полях эта площадь постепенно уменьшается до 2150 точек по мере роста заблаговременности; аналогичная тенденция просматривается и в максимальных значениях (снижение от 19812 до 18426). В этом можно увидеть тенденцию к росту положительной асимметрии распределения прогностических площадей: правая ветвь относительно максимальной ординаты постепенно вытягивается. Такие явления отражают основное качество большинства статистических моделей – сглаживание и соответствующее подавление экстремумов в значениях величин и в их площадных характеристиках с разбиением связных объектов на части.

Оценки качества модели обобщенного распределения Парето. Проверка статистической значимости оцениваемых параметров масштаба и формы должна предваряться оценкой качества моделирования. Следует иметь в виду, что некоторые метрики качества моделей зависят (явно или неявно) от параметров модели, вследствие чего зависят и от методов оценивания этих параметров. Так, информационный критерий Акаике рассчитывается с помощью функции правдоподобия, которая используется в численных методах оценок MLE и GMLE, поэтому является побочным продуктом оценок параметров и зависит от особенностей рассчитанной кросс-ковариационной матрицы этих параметров. Но поскольку на выходе методов в наличии оказываются оценки параметров (правда, некоторые методы иногда выдают либо неестественное число, либо константу отсутствия), то с их помощью можно проверить соответствие аналитической функции плотности гистограмме исходных наблюдений с помощью безотказного критерия Хи-квадрат.

Рассмотрим критерии Акаике и Хи-квадрат в приложении к выборкам экстремумов, характеристики которых отражены в табл. 1. В табл. 2 собраны оценки параметров масштаба и формы распределения Парето максимальных размеров областей сплошных осадков, наблюдаемых в теплый период 2017 г. в полях радиолокатора Смоленск и в полях наукастинга. В табл. 3 собраны значения в градациях гистограмм, построенных по данным табл. 2. К оценкам параметров добавляются границы ДИ с уровнем доверия 95 %.

Значения параметра масштаба ниже не комментируем, напомним лишь, что на графике плотности распределения Парето обратная величина масштаба есть значение плотности распределения в условном нуле, в данном случае в точке 625. Сосредоточимся на параметре формы. Видно, что распределение обладает тяжелым хвостом – все оценки параметра формы положительные вместе с доверительными интервалами (за малозначительным исключением в столбце наблюдений obs, где нижняя граница ДИ равна -0.043).

Таблица 2. Критерии качества моделирования общим распределением Парето максимальных площадей объектов и оценки параметров масштаба и формы с границами их доверительных интервалов в приложении к полям осадков по наблюдениям (obs) ДМРЛ-С Смоленск (RUDL) в теплый период года (май-сентябрь 2017 г.) и в прогностических полях (30-120 мин) по модели STEPS. Метод оценки параметров – обобщенное максимальное правдоподобие (GMLE)

Table 2. Quality criteria for modeling (Akaike and Chi-square) by the general Pareto distribution of the maximum areas of objects, and estimates of the scale and shape parameters with their confidence intervals according to observations (obs) of DMRL-S radar of Smolensk (RUDL) in the warm season (May - September 2017) and in STEPS nowcasting fields (30-120 min). The parameter estimation method is the generalized maximum likelihood (GMLE)

	obs 30 60 90 120
GML_Akaike	1661 2150 2225 2364 2352
GML_XI2	9.001 11.797 28.153 16.109 16.163
GML_scl_1	1471 1518 1531 1475 1564
GML_scl_2	2468 2270 2261 2155 2279
GML_scl_3	3464 3023 2991 2834 2995
GML_shp_1	-0.043 0.025 0.049 0.080 0.053
GML_shp_2	0.336 0.326 0.337 0.356 0.333
GML_shp_3	0.715 0.627 0.625 0.632 0.613

Примечание. Метод оценки параметров обобщенного распределения Парето –GMLE, Akaike – критерий Акаике, XI2 - критерий Хи-квадрат, scl - оценка параметра масштаба, shp – оценка параметра формы. Числа при оценках параметров: 1 – нижняя (2.5%-ная) граница ДИ, 2 – оценка параметра, 3 – верхняя граница (97.5%-ная) ДИ. Красным цветом выделены локальные максимумы в строках критериев качества.

Таблица 3. Количество случаев в градациях гистограммы распределения максимальных площадей по данным, соответствующим табл. 2

Table 3. The number of cases in the histogram bins of the distribution of maximum areas according to data similar to Table 2

Наблюдения/ прогнозы (STEPS-мин)	Градации гистограммы
	1 2 3 4 5 6 7 8 9 10 11
ДМРЛ-С RUDL	38 23 10 7 5 3 3 2 0 1
STEPS-030 STEPS-060 STEPS-090 STEPS-120	55 29 10 9 7 3 2 3 0 1 57 30 10 10 4 6 0 5 0 0 1 63 31 10 9 5 7 1 3 1 1 63 30 9 10 6 6 0 3 2 1

Примечание. Исходные данные: ДМРЛ-С Смоленск, теплый период, порог Парето равен 625. Эквидистантные градации выделены методом Стурджеса.

Качество моделирования обобщенного распределения Парето по критерию Акаике и по Хи-квадрат падает при переходе к полям прогноза 30 мин почти одинаково – в 1.3 раза по сравнению с качеством моделирования наблюдений. В обоих рядах значений критериев имеются "локальные" максимумы – слабый – на 90 мин (Акаике) и заметный на 60 мин (Хи-квадрат). Локализация и значение максимума Хи-квадрат не нарушают общего хода критериев. Трудно сказать, чем вызваны такие выбросы, так как на первый взгляд неясно, чем отличаются значения в градациях гистограммы на строке STEPS-060 в табл. 3 от значений в других строках.

Несмотря на то, что в методе GMLE параметр формы имеет априорное бета-распределение, ограничивающее значения параметра формы интервалом [-0.5, +0.5], апостериорное распределение может само иметь другие параметры, приводящие к оценкам доверительных интервалов, превышающих априорные границы параметра формы. Среди оценок параметра формы методом МСМС также встречаются значения, выходящие за указанные границы. В табл. 2 верхние границы доверительных интервалов заметно превышают 0.5, что оказывается вполне обычным явлением и для данных по остальным локаторам.

Оценка и ранжирование методов с помощью критерия Хи-квадрат. Сведем в отдельную таблицу строки характеристик для параметра формы, реорганизовав их для удобства сопоставления и добавив значения критерия Хи-квадрат (табл. 4). Степени свободы, оцененные по количеству градаций гистограмм и учитываемые при определении критических значений критерия Хи-квадрат, расположились для обоих периодов года в интервале между 6 и 12. В таблицах [1] 5%-ные критические значения для указанных в скобках степеней свободы равны 12.592 (6), 16.919 (9), 18.307 (10), 19.675 (11) и 21.026 (12). Значения критерия Хи-квадрат на столбце прогноза 60 мин имеют характер явного выброса по сравнению со значениями этого критерия в остальных столбцах. Обратим внимание на ранги методов оценок параметров. Пиковые объекты в полях наблюдений и прогнозов на всех сроках моделируются по критерию качества лучше всего методом Байеса, хуже всего методом L-моментов, при этом и оценки параметра формы оказываются наибольшими в основном по методу Байеса и наименьшими – по методу L-моментов. Все оценки параметра формы и доверительные интервалы, кроме нижних границ ДИ для полей наблюдений, являются положительными, что однозначно подтверждает "тяжелохвостость" распределения Парето тех областей сплошных осадков, площадь которых превышает 625 узлов сетки (напомним, что эквивалентно квадрату 50×50 км).

Сопоставив оценки критерия Хи-квадрат, можно утверждать, что все значения критерия, кроме прогнозов на 60 мин, не отклоняют гипотезы применимости модели распределения Парето к выделенным объектам.

Таблица 4. Критерий Хи-квадрат и оценки параметра формы с границами доверительных интервалов по данным, соответствующим табл. 2

Table 4. Chi-square test and shape parameter estimates with confidence interval boundaries based on data corresponding to Table 2

	obs 30 60 90 120
GML_XI2	9.001 11.797 28.153 16.109 16.163
MLЕ_XI2	10.21 13.354 30.346 16.923 17.769
Lmo_XI2	10.777 14.121 32.105 18.581 19.092
Bay_XI2	7.773 11.125 26.734 15.132 15.404
лучший худший	Bay_XI2 Bay_XI2 Bay_XI2 Bay_XI2 Bay_XI2 Lmo_XI2 Lmo_XI2 Lmo_XI2 Lmo_XI2 Lmo_XI2
GML_shp_1	-0.043 0.025 0.049 0.08 0.053
MLE_shp_1	-0.044 0.016 0.041 0.075 0.045
Lmo_shp_1	-0.007 0.034 0.047 0.069 0.052
Bay_shp_1	-0.008 0.037 0.063 0.11 0.091
больший меньший	Bay/Lmo Bay/Lmo Bay/GML Bay Bay MLE/GML MLE MLE LMO/MLE MLE
GML_shp_2	0.336 0.326 0.337 0.356 0.333
MLE_shp_2	0.294 0.289 0.309 0.34 0.304
Lmo_shp_2	0.278 0.274 0.29 0.312 0.284
Bay_shp_2	0.343 0.328 0.344 0.396 0.354
больший меньший	Bay/GML Bay/GML Bay/GML Bay Bay Lmo Lmo Lmo/MLE Lmo Lmo/MLE
GML_shp_3	0.715 0.627 0.625 0.632 0.613
MLE_shp_3	0.632 0.561 0.577 0.604 0.562
Lmo_shp_3	0.506 0.496 0.492 0.503 0.47
Bay_shp_3	0.831 0.681 0.676 0.731 0.698
больший меньший	Bay Bay Bay Bay Bay Lmo Lmo Lmo Lmo Lmo

Примечание. Исходные данные: теплый период 2017 г., порог Парето равен 625. Использованы все методы оценки параметров – GMLE, MLE, L-moment, Bayesian. Числа при оценках параметров: 1 – нижняя (2.5%-ная) граница ДИ, 2 – оценка параметра, 3 – верхняя граница (97.5%-ная) ДИ. Зеленым цветом выделены значения критерия Хи-квадрат, не противоречащие принятию гипотезы о правомерности модели обобщенного распределения Парето. Красным цветом выделены случаи превышения значений критерия Хи-квадрат для 95%-ного уровня значимости (локальные максимумы на прогнозе 60 мин). Косым шрифтом выделены ранги методов по принципу наибольшего и наименьшего значений соответствующих характеристик.

Ранжирование методик по критерию Хи-квадрат дает следующую информацию. Наилучшее качество моделирования для данного локатора, заданного периода года и при указанном пороге Парето обеспечивает байесовский метод, к наихудшему качеству моделирования приводит использование метода L-моментов. Заметим, что в других случаях и для других локаторов ранги могут быть иными.

Оценка с помощью критерия Акаике. Если сфокусироваться на одном параметре, скажем, параметре формы, то можно получить сравнительную таблицу для всех локаторов (табл. 5). Для краткости исключим доверительные интервалы и рассмотрим данные теплого периода с порогом Парето 625, используя критерий Акаике.

Таблица 5. Объемы выборок максимальных площадей осадков (peaks) и оценка качества моделирования по критерию Акаике обобщенного распределения Парето (при пороге 625) для радиолокаторов Центрального федерального округа по данным наблюдений (obs) и прогнозирования (30-120 мин) в теплый период 2017 г. Метод оценки параметров – GMLE

Table 5. Sample volumes of maximum precipitation areas (peaks) and evaluation of the quality of modeling of the generalized Pareto distribution (GPD) (625 grid points threshold) by the Akaike criterion for radars of the Central Federal District of Russia, observations (obs) and forecasting (30-120 min) in the warm period of 2017. Method of the parameter estimates is GMLE

ДМРЛ-С		Наблюдения ДМРЛ-С	Срок прогноза, модель STEPS, мин
		obs	30	60	90	120
Курск, RAKU	peaks Akaike	86 117 118 121 123 1496 2098 2123 2179 2219
Тула, RATL	peaks Akaike	79 117 117 128 132 1463 2103 2115 2298 2367
Воейково, RAVO	peaks Akaike	89 133 141 145 157 1620 2381 2534 2609 2816
Брянск, RUDB	peaks Akaike	76 90 96 99 99 1360 1609 1709 1760 1768
Курск, RUDK	peaks Akaike	96 130 142 145 155 1740 2346 2564 2631 2800
Смоленск, RUDL	peaks Akaike	92 119 123 131 130 1661 2150 2225 2364 2352
Н.Новгород, RUDN	peaks Akaike	84 122 129 141 147 1713 2415 2379 2588 2689
Валдай, RUWJ	peaks Akaike	85 119 128 138 152 1582 2202 2359 2533 2772

Как было отмечено ранее [5], количество объектов, в том числе и количество объектов значительной площади в прогностических полях, заметно вырастает уже на первых шагах прогнозирования (до 30 мин). Это количество растет и на остальных сроках, хотя и не столь заметно. Рассмотрим, к примеру, отношения количества объектов на шаге прогноза 120 мин к количеству объектов в наблюдениях (123/86 = 1.43 и т. д) и сопоставим с аналогичным отношением значений критерия Акаике (2219/1496 = 1.48 и т. д.). Обнаруживается немного удивительный факт: чем больше выборка, тем примерно в той же пропорции падает качество моделирования этой выборки распределением Парето. Так, для радиолокатора Брянск (RUDB) наблюдается почти точное совпадение этих пропорций, составляющих рост объема выборки и падение качества моделирования в 1.3 раза. Значения критерия Акаике для данных локатора RUDB оказываются наилучшими по всем столбцам таблицы.

Оценка параметра формы. По аналогичной схеме составим таблицу значений параметра формы, оцененного всеми четырьмя методами (табл. 6).

Таблица 6. Оценки параметра формы обобщенного распределения Парето (при пороге 625) для радиолокаторов Центрального федерального округа по данным наблюдений (obs) и прогнозирования (30-120 мин) в теплый период 2017 г.

Table 6. Estimates of the shape parameter of the generalized Pareto distribution (at a threshold of 625) for the Central Federal District radars based on observational data (obs) and forecasting (30-120 min) in the warm period of 2017

ДМРЛ-С		Наблюдения ДМРЛ-С	Срок прогноза, модель STEPS, мин
		obs	30	60	90	120
Курск, RAKU	GMLE MLE Lmom Bayes	0.428 0.405 0.413 0.390 0.366 0.472 0.419 0.434 0.394 0.353 0.379 0.354 0.360 0.336 0.312 0.681 0.560 0.628 0.533 0.447
Тула, RATL	GMLE MLE Lmom Bayes	0.278 0.360 0.344 0.376 0.371 0.161 0.342 0.316 0.370 0.363 0.181 0.307 0.290 0.329 0.324 0.179 0.424 0.361 0.432 0.436
Воейково, RAVO	GMLE MLE Lmom Bayes	0.228 0.294 0.307 0.326 0.314 0.074 0.241 0.264 0.294 0.280 0.094 0.238 0.255 0.275 0.267 0.109 0.260 0.273 0.341 0.307
Брянск, RUDB	GMLE MLE Lmom Bayes	0.327 0.303 0.328 0.348 0.322 0.271 0.236 0.284 0.319 0.276 0.266 0.237 0.275 0.301 0.268 0.342 0.280 0.355 0.397 0.311
Кострома, RUDK	GMLE MLE Lmom Bayes	0.284 0.331 0.328 0.303 0.361 0.199 0.299 0.298 0.259 0.350 0.206 0.281 0.281 0.252 0.316 0.234 0.355 0.336 0.288 0.417
Смоленск, RUDL	GMLE MLE Lmom Bayes	0.336 0.326 0.337 0.356 0.333 0.294 0.289 0.309 0.340 0.304 0.278 0.274 0.290 0.312 0.284 0.357 0.322 0.354 0.385 0.336
Н.Новгород, RUDN	GMLE MLE Lmom Bayes	0.000 0.000 0.253 0.308 0.322 0.000 -0.002 0.165 0.262 0.287 0.038 0.153 0.180 0.250 0.267 0.038 0.135 0.167 0.307 0.320
Валдай, RUWJ	GMLE MLE Lmom Bayes	0.287 0.249 0.262 0.291 0.327 0.190 0.153 0.184 0.237 0.296 0.202 0.169 0.194 0.234 0.276 0.203 0.162 0.196 0.247 0.318

Примечание. Параметры распределения Парето оценены с помощью обобщенного метода максимального правдоподобия (GMLE), максимального правдоподобия (MLE), L-моментов (Lmom) и моделирования цепей Маркова с помощью алгоритма Метрополиса-Гастингса (Bayesian).

Здесь повсюду наблюдается "тяжелохвостость" распределения, кроме области обзора локатора Нижний Новгород, где значения параметра формы близки к нулю и слабо отрицательны для 30 мин. При этом оценки Акаике в табл. 5 на строке RUDN никак не выделяются на общем фоне аналогичных оценок для остальных локаторов.

При переходе к прогнозу на 30 мин величина параметра формы может как увеличиваться (RATL, RAVO, RUDK), так и уменьшаться (RAKU, RUDB, RUDL, RUWJ), причем для всех методов оценки параметров. При переходе по прогностическим срокам возможны оба варианта: и уменьшения, и увеличения.

Сравнение оценок параметра формы, полученных различными методами, можно провести по минимальным и максимальным значениям в соответствующих столбцах. Так, например, в данных для локатора RAKU в столбце наблюдений максимальное значение параметра формы равно 0.681 и получено методом Байеса, а минимальное значение равно 0.379 и получено методом L-моментов.

Ранжирование методов по всем столбцам таблицы, с учетом близких значений (когда ранги приравниваются, а методы отмечаются совместно), привело к следующей статистике.

Максимальные значения параметра формы поставляются методами Байеса (21 случай) и обобщенным методом максимального правдоподобия (19 случаев). Минимальные значения параметра формы получены в основном методом L-моментов (30 случаев) и лишь в 10 случаях – методом максимального правдоподобия.

Причины сходства методов, порождающих максимальные оценки параметра формы, вполне очевидны: оба основаны на байесовом подходе либо целиком (Bayesian), либо в "усеченном" виде (GMLE). Как правило, оценки методом MLE ниже оценок родственного метода GMLE. Имеются случаи максимального расхождения оценок MLE и GMLE: в данных наблюдений RATL (0.161 и 0.278), RAVO (0.074 и 0.228), RUDK (0.199 и 0.284), RUWJ (0.190 и 0.287), а также в двух-трех столбцах прогностических данных.

Выходные графические данные. Графическая продукция раскрывает другие стороны моделирования значительных площадей осадков с помощью распределения Парето. Графика для каждого локатора и каждого порога Парето сохраняется в одном файле на 20 страницах (слайдах) в следующем составе. Для каждого из четырех методов зарезервировано по пять страниц, соответствующих наблюдению и четырем прогнозам. На каждой из страниц для методов GMLE, MLE и Lmoments содержатся четыре панели:

(1) квантильная диаграмма "модель – наблюдение";

(2) квантильная диаграмма с линией регрессии "наблюдение – имитация";

(3) плотность распределения наблюдения (в ядерном сглаживании) и модельная плотность Парето;

(4) период повторяемости с доверительными интервалами.

На каждой из страниц для метода Bayesian содержатся четыре панели:

(1) и (2) – апостериорная плотность распределения параметра масштаба и формы соответственно (в ядерном сглаживании);

(3) и (4) – "следы" реализаций марковской последовательности для апостериорных плотностей распределения параметров масштаба и формы (10 тысяч шагов, с указанием уровня "выгорания" в burn in = 500 шагов).

Для оценки качества моделирования накапливались и частично анализировались все критерии, содержащиеся в модуле fevd (AIC, BIC и DIC), дополнительно рассчитывался критерий Хи-квадрат. Как указывалось выше, ни один из встроенных критериев качества не применяется ко всем четырем методам оценки параметров. В пользу критерия Хи-квадрат можно добавить следующее. Информационный критерий Акаике рассчитывается по формуле 4k – 2ln(L), где L – максимальное значение функции правдоподобия; k – количество параметров в модели (в нашем случае всегда k=2). Но в предположении нормального распределения ошибок модели критерий Акаике совпадает с точностью до аддитивной константы с числом c²+2k.

Оценка качества моделирования графическими средствами. Проиллюстрируем табличные данные для локатора Смоленск (RUDL) (табл. 1–4) примерами выходной графической продукции модуля fevd(). Ограничимся наблюдениями и прогнозами на 90 минут.

На рис. 1 и 3 демонстрируются квантильные диаграммы методов GMLE, MLE и L-moments, на рис. 2 и 4 – апостериорные плотности распределения параметров масштаба и формы вместе с их следами случайного блуждания по методу Bayesian.

Квантильные диаграммы строятся двумя способами:

1) между наблюдениями и теми модельными данными, которые восстанавливаются по соответствующим квантилям обратным преобразованием распределения Парето с оцененными параметрами масштаба и формы;

2) между наблюдениями и имитационными данными, сгенерированными датчиком случайных чисел, настроенным на распределение Парето с указанными оценками параметров.

На рис. 1 и 3 первый набор диаграмм (obs-model) размещен в левом столбце панелей, а второй набор диаграмм (simul-obs) – в правом столбце панелей.

Первое и, пожалуй, основное замечание следует сразу сделать об объектах самых больших размеров, которые были сохранены в выборках "пиков над порогом", хотя иногда такие объекты представляли собой области осадков, выходящие за пределы кругов обзоров локаторов, нарушая физическую, а отсюда и статистическую корректность выборки.

obs-model qq-diagram	simul-obs qq-diagram
GMLE, XI2 = 9.00
MLE, XI2 = 10.21
L-moments, XI2 = 10.78

Рис. 1. Квантильные диаграммы для сопоставления данных наблюдений и данных моделирования распределением Парето с оцененными параметрами масштаба и формы (obs-model, левый столбец панелей), а также имитационных данных со стохастическим генерированием выборок на оцененных параметрах распределения Парето (simul-obs, правый столбец панелей). Порог распределения Парето - 625 точек. Методы оценки параметров: обобщенное максимальное правдоподобие (GMLE), максимальное правдоподобие (GMLE), L-моменты (L-moments). Данные радиолокатора Смоленск (RUDL), теплый период 2017 г.

Fig. 1. Quantile charts for comparing observational data and Pareto distribution simulation data with estimated scale and shape parameters (obs - model, left panel column), as well as simulation data with stochastic sampling on estimated Pareto distribution parameters (simul - obs, right panel column). Pareto distribution threshold - 625 points. Parameter estimation methods: generalized maximum likelihood (GMLE), maximum likelihood (GMLE), L - moments (L - moments ). Smolensk radar data (RUDL), warm period 2017.

Исключение такого рода объектов (цензурирование выборки) предусмотрено опцией max.size в функции FeatureFinder() библиотеки SpatialVx, с помощью которой были сформированы пространственные объекты для верификации [5], однако сделать это в автоматизированном режиме затруднительно. Необходимость цензурирования напрашивается при анализе квантильных графиков левого столбца рис. 1 и 3. Например, большее совпадение модельных и наблюдаемых данных могло быть достигнуто при исключении объектов размеров более 20 тысяч точек, что, кстати, эквивалентно квадрату со стороной, равной радиусу круга обзора локатора (250 км).

Рис. 2. Апостериорные плотности распределения параметра масштаба (левая верхняя панель) и параметра формы (правая верхняя панель), и следы случайного блуждания параметра масштаба (левая нижняя панель) и параметра формы (правая нижняя панель), построенные по байесовскому методу (Bayesian) с генерированием марковских цепей алгоритмом Метрополиса-Гастингса (Монте-Карло). Критерий Хи-квадрат = 8.3. Данные радиолокатора Смоленск (RUDL), теплый период 2017 г.

Fig. 2. Posterior distribution densities of the scale parameter (upper left panel) and shape parameter (upper right panel) and random walk traces of the scale parameter (lower left panel) and shape parameter (lower right panel) plotted by Bayesian method using Markov chain generation on the basis of the Metropolis-Hastings (Monte Carlo) algorithm. Chi-squared test = 8.3. Smolensk radar data (RUDL), warm period 2017.

obs-model qq-diagram	simul-obs qq-diagram
GMLE, XI2 = 16.109
MLE, XI2 = 16.923
L-moments, XI2 = 18.581

Рис. 3. То же, что на рис. 2, кроме наклонного жирного текста, который следует заменить на "Прогнозы модели STEPS на 90 мин, зона обзора и начальные данные радиолокатора Смоленск (RUDL), теплый период 2017 г."

Fig. 3. The same as in fig. 2, except for italic bold text, which should be replaced with "STEPS model forecasts for 90 min, field of view and initial data of the Smolensk (RUDL) radar, warm period 2017."

Несмотря на наличие выбросов в квантильных графиках левого столбца, удается получить, хотя и с широким веером доверительных интервалов, приемлемое совпадение на имитационных данных. Например, на рис. 1 заметное отклонение от линии совпадения одного-двух экстремумов на диаграммах левого столбца компенсируется вполне удачным имитированием на диаграммах правого столбца для методов GMLE и MLE. Этого, к сожалению, нельзя сказать о данных для полей прогноза на 90 мин (рис. 3), качество моделирования которых по критерию Хи-квадрат заметно уступает качеству моделирования объектов наблюдения (рис. 1). В любом случае при цензурировании размеров объектов верхней границей (около 15–20 тысяч точек) диаграммы обеих столбцов окажутся более привлекательными.

Рис. 4. То же, что на рис. 2, кроме наклонного жирного текста, который следует заменить на "Прогнозы модели STEPS на 90 мин, зона обзора и начальные данные радиолокатора Смоленск (RUDL), теплый период 2017 г."

Fig. 4. The same as in figure 2, except for italic bold text, which should be replaced with "STEPS model forecasts for 90 min, field of view and initial data of the Smolensk (RUDL) radar, warm period 2017."

Апостериорные плотности распределения и следы случайного блуждания оценок параметров демонстрируются на рис. 2. Выведены данные для объектов в полях наблюдений и в полях прогнозов на 90 мин. На графиках следов выделена линия отжига, исключающая из построения апостериорных плотностей все реализации марковской цепи левее точки 500. Сравнивая рис. 2 и 4, отметим следующее. Моды плотностей для масштаба и для формы в прогностических выборках отклоняются в разные стороны – масштаб снижется примерно на 500 точек (кривая плотности распределения Парето в "нуле" поднимается), а параметр формы увеличивается примерно с 0.4 до 0.5 (хвост распределения Парето "тяжелеет"). В частности, такое поведение кривой плотности распределения Парето предполагает, что при прогнозе на 90 мин модель STEPS в окрестности условного "нуля" (эквивалента порогу 625 точек) систематически завышает, далее занижает, а начиная с некоторого размера объектов начинает систематически завышать площади объектов. Рост параметра формы до 0.5 означает, в частности, что у распределения Парето для прогнозируемых объектов может остаться только математическое ожидание, а оценка дисперсии по таким данным потеряет смысл.

2.2. Выбор метода оценивания параметров и порога распределения Парето

Основываясь на полученных оценках формы и их доверительных интервалов, можно утверждать, что все четыре метода привели к согласованным выводам о параметре формы для порога 625.

Стандартным методом оценок является, несомненно, метод максимального правдоподобия (MLE), связанный, кстати говоря, с видоизмененным методом минимума Хи-квадрата [3]. Однако на малых выборках MLE может приводить к неестественным оценкам параметров, что было причиной предложения поправки в виде урезанного байесовского варианта (GMLE) [12].

Методы L-моментов привлекательны ввиду простоты расчетов и статистической робастности оценок, однако в [12] с помощью статистических экспериментов с большими квантилями было показано преимущество метода GMLE и перед методом L-moments в условиях выборок среднего объема и положительного параметра формы.

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

Ограничение одним методом GMLE лишь снижает исходную многопараметричность и неопределенность решаемой задачи, но оставляет широкий простор для предположений и допущений. Хотя оценка качества модели ограничена одним критерием Хи-квадрат, а доверительные интервалы – одним уровнем (95 %), следует иметь в виду разные методы расчетов доверительных интервалов для разных методов оценки параметров распределения Парето. Поэтому, не претендуя на предложение универсального интегрального критерий качества, мы делаем выводы исключительно в сравнительных терминах (больше-меньше) относительно произвольно выбранных пороговых значений.

На рис. 5 представлены гистограммы распределения размеров объектов, не меньших 625, 900, 1225 и 1600 точек, с нанесенными на гистограммы значениями плотности распределения Парето, соединенными линейными отрезками. В названиях панелей указаны объемы выборок, автоматически рассчитанное количество градаций, экстремальные значения и медиана размеров, а также искомые оценки параметров. Для удобства визуального сравнения пары гистограмм одного ряда (наблюдение – прогноз) приводятся к единой градуировке с помощью предварительных построений "заглушенных" гистограмм для объединения двух выборок. В функции hist() количество k градаций определяется с помощью правила Стурджеса: k = 1 + [lg₂(n)], где n – объем выборки.

Хорошо видны отличия и сходства между гистограммами и аппроксимирующими кривыми плотности распределения Парето. Как уже известно, в среднем количестве в полях прогнозов генерируется больше объектов, чем в полях наблюдений. Так как больший порог Парето выделяет подмножество из тех максимумов, которые отобраны и для меньшего порога, то аппроксимация выделенного подмножества распределением Парето необходимо увеличивает масштаб (уход вправо по хвосту) и меняет форму.

Выпишем оценки критерия Хи-квадрат для рис. 5:

625: χ² (RAKU) = 13.190, χ² (STEPS-60) = 14.721,

900: χ² (RAKU) = 13.106, χ² (STEPS-60) = 14.201,

1225: χ² (RAKU) = 13.420, χ² (STEPS-60) = 15.788,

1600: χ² (RAKU) = 14.099, χ² (STEPS-60) = 24.245.

Выпишем оценки критерия Хи-квадрат для рис. 6:

625: χ² (RAKU) = 10.849, χ² (STEPS-60) = 13.646,

900: χ² (RAKU) = 12.427, χ² (STEPS-60) = 10.810,

1225: χ² (RAKU) = 16.501, χ² (STEPS-60) = 10.982,

1600: χ² (RAKU) = 22.069, χ² (STEPS-60) = 15.124.

Визуальный анализ рис. 5 и 6 позволяет сделать следующий качественный вывод для моделирования данных из зоны обзора локатора Курск в период май–сентябрь 2017 г. Наиболее подходящим является пороговый размер объекта порядка 600–900 точек. Объекты размера порядка 1200 и выше моделируются явно неудовлетворительно. Аналогичный анализ для остальных локаторов дает основания распространить сделанный вывод на все использованные локаторы испытаний по территории ЦФО.

Подтверждение предыдущих выводов представим табл. 7, в которой собраны оценки Хи-квадрат для соответствующих данных со всех локаторов (столбцы RADAR) в сравнении с прогностическими данными заблаговременностью 60 мин (столбцы STEPS-60).

Наблюдение: локатор RAKU	Прогноз на 60 мин: STEPS
size ≥ 625 vol 87 - 119
size ≥ 900 vol 69 - 97
size ≥ 1225 vol 53 - 74
size ≥ 1600 vol 46 - 65

Рис. 5. Теплый период. Гистограммы и приближение распределением Парето размеров объектов в полях осадков в наблюдениях локатора Курск (RAKU, левый столбец) и в прогнозах (STEPS-60, правый столбец) на 60 мин. Порог Парето равен (сверху вниз) 625, 900, 1225 и 1600 точек поля. Значения размеров приведены в масштабе size/1000.

Fig. 5. Warm period. Histograms and Pareto distribution approximation of object sizes in precipitation fields in Kursk radar observations (RAKU, left column) and forecasts (STEPS-60, right column) for 60 min. The Pareto threshold is (from top to bottom) 625, 900, 1225 and 1600 field points. Dimensions are given in size/1000 scale.

Наблюдение: локатор RAKU	Прогноз на 60 мин: STEPS
size ≥ 625 vol 59 - 106
size ≥ 900 vol 51 - 81
size ≥ 1225 vol 44 - 68
size ≥ 1600 vol 41 - 55

Рис. 6. Холодный период. Гистограммы и приближение распределением Парето размеров объектов в полях осадков в наблюдениях локатора Курск (RAKU, левый столбец) и в прогнозах (STEPS-60, правый столбец) на 60 мин. Порог Парето равен (сверху вниз) 625, 900, 1225 и 1600 точек поля. Значения размеров приведены в масштабе size/1000.

Fig. 6. Cold period. Histograms and Pareto distribution approximation of object sizes in precipitation fields in Kursk radar observations (RAKU, left column) and forecasts (STEPS-60, right column) for 60 min. The Pareto threshold is (from top to bottom) 625, 900, 1225 and 1600 field points. Dimensions are given in size/1000 scale.

Таблица 7. Значения критерия Хи-квадрат для оценки качества приближения гистограммы распределением Парето с оцененными параметрами масштаба и формы при порогах 625, 900, 1225 и 1600 точек

Table 7. Chi-square test values for assessing the quality of histogram approximation by the Pareto distribution with estimated scale and shape parameters at thresholds of 625, 900, 1225, and 1600 points

РАДАР (cases; ndeg)	RADAR	STEPS-60 MIN
	625 900 1225 1600	625 900 1225 1600
Теплый период
RAKU (87-46;11) RATL (80-57;11) RAVO (90-65;11) RUDB (77-52;11) RUDK (97-66;10) RUDL (93-61;12) RUDN (85-61;13) RUWJ (86-61;11)	13.190 13.106 13.420 14.009 4.097 5.674 8.708 12.476 6.921 8.417 13.148 19.733 6.690 9.289 13.859 21.479 6.862 8.432 13.818 18.803 4.711 6.310 9.959 13.958 4.301 5.301 9.213 12.080 4.059 4.499 9.204 12.638	14.721 14.201 15.788 24.245 11.395 12.140 10.565 17.086 4.637 7.423 11.288 21.318 7.695 10.024 15.256 21.446 5.947 8.044 13.332 20.265 14.335 16.637 22.859 30.741 7.474 10.088 14.666 22.267 5.220 7.727 13.841 21.949
Холодный период
RAKU (59-41;12) RATL (48-25;13) RAVO (54-29;10) RUDB (46-25;10) RUDK (41-18;10) RUDL (47-25;7) RUDN (41-20;11) RUWJ (27-16;12)	10.849 12.427 16.501 22.069 10.775 10.747 11.624 11.623 10.296 11.282 15.997 23.012 7.617 8.678 11.632 9.365 11.502 13.442 13.641 16.047 3.604 3.020 2.056 1.888 12.630 11.117 10.904 12.481 8.370 10.298 14.593 24.938	13.646 10.810 10.982 15.125 13.599 12.759 10.450 7.272 8.947 9.096 15.831 23.653 7.689 7.192 10.107 14.690 12.697 12.139 15.109 22.009 12.928 4.424 3.651 2.909 6.423 5.255 7.355 10.534 13.222 12.296 13.484 18.493

Примечание. В столбце РАДАР рядом с идентификатором в скобках указаны (объем выборки наблюдений для порога 625 точек и для порога 1600 точек; оценка числа степеней свободы). Желтый цвет – Хи-квадрат в критической области, красный цвет – рост качества модели с ростом порога Парето.

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

Красным цветом выделены значения в зоне обзора локатора Смоленск в холодный период, отражающие, на первый взгляд, одно из важнейших условий второй теоремы об экстремумах, пороговую устойчивость: чем больше порог, тем точнее данные моделируются распределением Парето. В реальных выборках ограниченного объема, стремительно уменьшающегося с ростом порога, данное явление следует признать редкой удачей. Рассмотрим подробнее соответствующие данные. На рис. 7 собраны гистограммы размеров максимальных объектов в полях наблюдений и прогнозов на 60 мин размером 625, 900, 125 и 1600 точек. Видно, что уменьшение с ростом порога значения критерия Хи-квадрат обеспечивается все более точным приближением кривой плотности распределения Парето первой градации гистограмм.

Наблюдение: локатор RUDL	Прогноз на 60 мин: STEPS
size ≥ 625 vol 47 - 127
size ≥ 900 vol 38 - 91
size ≥ 1225 vol 29 - 74
size ≥ 1600 vol 24 - 58

Рис. 7. Холодный период. Гистограммы и приближение распределением Парето (синяя линия) размеров объектов в полях осадков в наблюдениях локатора Смоленск (RUDL, левый столбец) и в прогнозах (STEPS-60, правый столбец) на 60 мин. Порог Парето равен (сверху вниз) 625, 900, 1225 и 1600 точек поля. Значения размеров приведены в масштабе size/1000.

Fig. 7. Cold period. Histograms and Pareto distribution approximation (blue line) of object sizes in precipitation fields in Smolensk locator observations (RUDL, left column) and in forecasts (STEPS-60, right column) for 60 min. The Pareto threshold is (from top to bottom) 625, 900, 1225 and 1600 field points. Dimensions are given in size/1000 scale.

Полный набор ситуаций, породивших искомые пики над порогами, представлен в табл. 8 для порога 625 точек; при росте порога выборка составляется из размеров этих же объектов.

Таблица 8. Пространственные, статистические и временные характеристики ситуаций с объектами размером не менее 625 точек в сетке двухкилометрового разрешения в зоне радиолокатора Смоленск. Холодный период ноябрь 2017 – март 2018 г. Объекты выделены изолинией 1 мм/ч после пространственного осреднения с радиусом 9 точек

Table 8. Spatial, statistical and temporal characteristics of situations with objects with a size of at least 625 points in a two-kilometer resolution grid observed in the Smolensk radar zone during the cold period November 2017 - March 2018 The objects are marked with a 1 mm/h isoline after spatial averaging with a radius of 9 points

Примечание. Номера ситуаций указаны в столбце sit. Далее следуют квартильные характеристики размеров объектов в данной ситуации(min, q25, med, q75 и max). Количество выделенных объектов обозначено именем valid, длительность ситуации в минутах – mins; срок первого поля ситуации –d_time_start; дата и срок последнего поля ситуации - d_time_stop. Если не указаны часы и минуты в столбцах дат и времени, то это соответствует времени суток 00:00.

В табл. 8 под номерами 32 и 33 показаны две ситуации, в которых максимальный объект оказался единственным. При этом под номером 33 ситуация имела продолжительность 60 мин, а объект имел размер 13044 точки, что эквивалентно квадрату со стороной около 230 км (почти четверть от зоны обзора локатора). Однако самый большой по размеру объект возник в ситуации 18, продолжавшейся 1280 мин (почти сутки). При этом в данной ситуации в 128 полях радиолокационных наблюдений имелось в совокупности 125 объектов размером не менее 632 точки. Вполне возможно, что в реальной синоптической ситуации общее количество объектов, сконструированных пространственным осреднением в 9 узлов сетки и изогиетой 1 мм, представляет одну систему осадков фронтального типа.

Отмеченное выше благоприятное свойство набора экстремальных величин в столбце max обусловлено как неизменным количеством эквидистантных градаций (всего 6 интервалов) для всех порогов, так и плотным расположением сорока двух первых значений вариационного ряда в диапазоне от 632 до 5617 (рис. 8). Если первое обстоятельство искусственно, то второе обстоятельство реально отражает пороговую устойчивость анализируемой выборки.

Рис. 8. График вариационного ряда, составленного из набора сорока семи экстремумов табл.8 (столбец max).

Fig. 8. Graph of ordered statistics, made up of a set of forty-seven extrema in Table 8 (column max).

2.3. Интегральные оценки качества с помощью параметров распределения

Продолжим сопоставление параметров распределения Парето для максимальных объектов в ситуациях наблюдений и прогнозов на основе данных радиолокатора Курск (RAKU) (табл. 9).

Таблица 9. Количество ситуаций (#situat), оценки параметров масштаба (scale), формы (shape) распределения Парето размеров объектов с порогами 625 и 900 точек, доля пересечения доверительных интервалов оценок параметров (intersect). Объекты извлекаются из ситуаций в последовательностях полей наблюдений (RAKU) и прогнозов (STEPS) с заблаговременностью 30, 60, 90 и 120 мин. Данные наблюдений и прогнозов относятся к двум периодам: теплому (май-сентябрь 2017 г.) и холодному (ноябрь 2017 – март 2018)

Table 9. Number of situations (# situations), estimates of scale parameters (scale), shape (shape) of the Pareto distribution of object sizes with thresholds of 625 and 900 grid points, the intersection ratio (intersect) of confidence intervals of parameter estimates. Objects are extracted from situations in series of fields of observations (RAKU) and forecasts (STEPS) with lead times of 30, 60, 90 and 120 min. Observation and forecast data refer to two periods: warm - May-September 2017 and cold - November 2017 - March 2018

Порог		Теплый период			Холодный период
		scale	shape	#situat	scale	shape	#situat
Заблаговременность 30 мин
625	RAKU	1956 0.428 86			2617 0.294 59
	STEPS	1859 0.420 117			2105 0.410 87
	intersect	0.80 0.84			0.75 0.78
900	RAKU	2497 0.362 69			3415 0.000 51
	STEPS	2504 0.321 94			2833 0.324 70
	intersect	0.74 0.74			0.68 0.47
Заблаговременность 60 мин
625	STEPS	1979 0.413 118			1875 0.457 106
	intersect	0.83 0.85			0.50 0.68
900	STEPS	2621 0.320 97			2891 0.320 81
	intersect	0.77 0.76			0.63 0.44
Заблаговременность 90 мин
625	STEPS	1958 0.409 121			1342 0.558 129
	intersect	0.83 0.83			0.23 0.54
900	STEPS	2268 0.368 103			2248 0.395 96
	intersect	0.73 0.76			0.38 0.33
Заблаговременность 120 мин
625	STEPS	2062 0.382 123			1338 0.548 142
	intersect	0.79 0.80			0.21 0.54
900	STEPS	2250 0.364 107			1837 0.459 112
	intersect	0.68 0.72			0.20 0.23

Рассмотрим результат прогноза на 60 мин. Добавим оценки ошибок, в таблицы не включенные. Из наблюдений в теплый период выделено 86 ситуаций; оценка параметра масштаба и стандартная ошибка равны 1956 и 472.9, соответственно; оценка параметра формы и стандартная ошибка равны 0.428 и 0.212. Напомним, границы доверительного интервала в методе GMLE определяются как "оценка ± 1.96*ошибка". Аналогичные числа для прогноза на 60 мин таковы: выделено 118 ситуаций; оценка параметра масштаба и стандартная ошибка равны 1979 и 391.1 соответственно; оценка параметра формы и стандартная ошибка равны 0.413 и 0.180.

Так как интегральная оценка качества построена на доверительных интервалах, обсудим их более подробно. Запишем два сравниваемых ДИ с помощью границ в виде (L₁, U₁) и (L₂, U₂) соответственно. Доля пересечения (intersect), визуально очевидная, рассчитывается по формуле с условием: intersect = (min(U₁,U₂) – max(L₁,L₂))/(max(U₁,U₂)-min(L₁,L₂)); при intersect < 0 присваивается нуль. Для параметра масштаба в выбранном случае intersect = 0.83, для параметра формы intersect = 0.85. Используя оценки ошибок, убедимся, что доверительные интервалы параметров прогностических полей являются собственными подмножествами соответствующих интервалов для параметров наблюдений, а это значит, что надежность оценки параметров выше в моделях для прогностических наборов. В холодный период 2017–2018 гг. были получены следующие результаты: доля пересечения для параметра масштаба равно 0.50, а для параметра формы – 0.68.

Статистическая значимость оценки intersect не оценивалась, выводы делаются в виде простого порядкового сравнения (больше-меньше) c произвольным представлением об уровне информативности и полезности. При этом информация о границах ДИ, утраченная в доле пересечения, может быть существенной при интерпретации результатов верификации. Действительно, при фиксированной доле пересечения положения ДИ относительно друг друга могут быть двух видов. Во-первых, ДИ параметра прогностического распределения может частично выходить левее или правее ДИ того же параметра для набора наблюдений. В данной ситуации теряется информация о возможных систематических особенностях прогностической системы. Во-вторых, один интервал может быть строго включен в другой интервал. Не зная принадлежности вложенного ДИ к прогнозу или к наблюдению, мы не сможем оценить, какой из двух наборов данных обеспечил более надежную оценку параметра.

В дополнение к данным табл. 9 приведем оценки параметра формы для площадей 1225 и 1600, добавив стандартную ошибку оценки параметра для прогноза на 60 мин и для соответствующих наблюдений. Зафиксированная ниже перемена знака параметра формы примерно на пороге Парето в 1225 точек оказывается характерной для большинства использованных локаторов. При изменении порога Парето в последовательности 625, 900, 1225 и 1600 объемы выборок из наблюдений менялись в последовательности 86, 69, 53 и 45. Объемы прогностических выборок, соответственно, равнялись 117, 94, 72 и 63. Рост объема прогностических выборок примерно одинаков для всех порогов: 117/86  » 1.36, 1.38, 1.36 и 1.4. Отметим, что такое постоянство не наблюдается во всех случаях.

Оценки параметра формы для наблюдений равнялись 0.428, 0.362, 0.000, -0.126; для прогнозов – 0.413, 0.320, -0.041, -0.172. Параметр формы для обоих наборов меняет знак на пороге в 1225 точек.

Последовательности стандартных ошибок выглядят следующим образом. Для наблюдений – 0.220, 0.295, 0.110 и 0.179; для прогнозов – 0.183, 0.219, 0.133, 0.109. Ошибки на переходах от порога 625 к 900 растут в обоих случаях, на переходах к порогам 1225 и 1600 такой согласованности не наблюдается. Согласно общей концепции, обсужденной в [5], при росте порога качество моделирования экстремумов должно расти благодаря асимптотическому приближению, но при этом из-за редеющих выборок ошибки оценок параметров должны также расти. На примере локатора Курск поведение ошибок укладывается в теоретические предположения только на переходе от первого порога ко второму. Критерий Хи-квадрат на этих переходах уменьшается (хотя и незначительно) в полном соответствии с теорией (табл. 7): для наблюдений от 13.2 до 13.1, для прогнозов от 14.7 до 14.2. Влияние роста объема выборки на ошибку оценки параметра можно усмотреть в том, что ошибки оценок параметра формы для прогнозов меньше соответствующих ошибок для наблюдений. В этом случае и доверительные интервалы для прогностических наборов более узкие по сравнению с доверительными интервалами для выборок наблюдений.

Далее, для идеального (в климатологическом смысле) прогноза и при удовлетворительном качестве модели распределения показатель intersect для каждого параметра распределения должен быть близким к единице. Но обратное не верно в общем: из равенства единице доли пересечения соответствующих доверительных интервалов совсем не следует идеальность прогноза в климатологическом смысле. Во-первых, одинаковые ранговые ряды любой исследуемой величины исключают из анализа время. Но здесь имеется надежная гарантия: максимумы не могут появиться в прогнозах далеко по времени от максимумов в наблюдениях по самой схеме наукастинга, которая стартует каждые десять минут и разнести далеко отстоящие экстремальные площади не способна принципиально. Во-вторых, доверительные интервалы могут представлять собой ворота, вылезающие за пределы графических рамок. Здесь гарантия состоит в такой узости доверительных интервалов, которую можно считать разумной и допустимой. Однако в случае тяжелых хвостов полезную информацию можно извлечь из расположения ДИ относительно нуля: качество прогноза тем выше, чем дальше относительно нуля смещены в положительную сторону оба ДИ – и для набора прогнозов, и для набора наблюдений. В этом случае размах ДИ менее значителен для выводов.

При всех возможных недостатках предложенного интегрального показателя качества, можно утверждать, что по данным табл.9 качество прогнозирования объектов существенного размера сроком на 60 мин в холодный период (2017–2018 гг.) в зоне обзора локатора Курск уступает качеству прогнозирования по данным того же локатора в теплый период.

Продолжим обобщение результатов верификации для всех локаторов, использованных в испытаниях 2017–2018 гг. Выделяем пороги Парето в 625 и 900 точек как наиболее информативные и дающие надежную основу для анализа данных. Для более устойчивой оценки параметров Парето по выделенным ситуациям в расчет допускались лишь выборки, в которых количество ситуаций было не меньше 20. Для обоих периодов такому условию в значительной степени не удовлетворяли локаторы Внуково (RAVN, зимой) и Валдай (RUWJ) – эти данные были полностью исключены.

Рассмотрим долю пересечения доверительных интервалов для оценок параметров масштаба и формы (табл. 10). Выделим некоторый уровень «неуспеха», например для значений intersect < 50%, и обозначим соответствующие ячейки красным цветом. В поведении доли пересечения имеются разнообразные особенности в зависимости от заблаговременности и порога Парето, частью систематические, частью случайные выборочные, здесь укажем лишь самые заметные из них, основываясь на значениях, выделенных красным цветом. Один из общих выводов заключается в том, что по полученным оценкам известное в прогностической практике более высокое качество прогноза осадков в холодный период не распространяется на прогнозы экстремальных по площади областей в зонах обзора всех локаторов.

Рассмотрим в самых общих чертах параметр масштаба. Оценка параметра влияет на значение плотности распределения в окрестности относительного нуля, и тем самым на вероятностные свойства объектов, имеющих площади, близкие к порогу Парето. Уменьшение масштаба увеличивает количество таких объектов, напротив, увеличение параметра приводит к уменьшению такого рода объектов в окрестности порога. Надо иметь в виду, что планомерное увеличение порога приводит к столь же планомерному увеличению параметра масштаба обобщенного распределения Парето.

Согласно расположению красных ячеек в столбцах масштаба прогнозирование по параметру более качественно в холодный период по данным локаторов RATL и RAVO, а в теплый период – по данным локаторов RAKU, RUDB и RUDL.

Обратимся к доле пересечения ДИ для параметра формы (shape, табл. 10). Сопроводим значения знаками (+/-) в зависимости от знака формы для пары (наблюдение / прогноз), знаком 0 – для параметра формы в интервале [-0.1, 0.1].

От прогностической системы желательно получить в первую очередь оценку способности сохранить знак параметра формы. По данным таблицы отчетливо видна зависимость знака формы от порога Парето. Так, в теплый период параметр при переходе от порога 625 к порогу 900 в общем снижается. Как отмечалось выше, при переходе к "забракованным" порогам 1225 и 1600 оценка параметра формы уходит к нулевому значению и даже к отрицательному. Эта граница между порогами дает дополнительные основания предполагать, что области с размером, превышающим 625 точек (примерно 50×50 км), наиболее подходят для анализа паретовости ("тяжелохвостости") как в наблюдениях, так и в прогнозах при данных объемах выборок.

Таблица 10. Доля пересечения доверительных интервалов (intersect) оценки параметра масштаба и параметра формы для порогов Парето 625, 900, 1225 и 1600 точек, для заблаговременностей 30, 60, 90 и 120 мин, для теплого и холодного периодов 2017–2018 гг. Красным цветом отмечены значения меньше 50 %. Пары символов (+/-/0), первый для наблюдения, второй – для прогноза, означают знаки параметра формы и попадание в интервал [-0.1, 0.1]

Table 10. Intersection ratios of confidence intervals (intersect) of the estimate of the scale parameter and the shape parameter for Pareto thresholds of 625, 900, 1225, and 1600 points, for lead times of 30, 60, 90, and 120 min, for the warm and cold periods of 2017–2018. Values less than 50% are marked in red. Pairs of symbols (+/-/0), the first for observation the second for forecasts, mark the parameter sign and the parameter falling into the range [-0.1, 0.1].

ДМРЛ	Порог/ Заблаговременность	Доля пересечения доверительных интервалов, %
		МАСШТАБ				ФОРМА
		Теплый период		Холодный период		Теплый период		Холодный период
		625	900	625	900	625	900	625	900
RAKU	30	80	74	75	68	84 (+ +)	74 (+ +)	78 (+ +)	47 (0 +)
	60	83	77	50	63	85 (+ +)	76 (+ +)	68 (+ +)	44 (0 +)
	90	83	73	23	38	83 (+ +)	76 (+ +)	54 (+ +)	33 (0 +)
	120	79	68	21	20	80 (+ +)	72 (+ +)	54 (+ +)	23 (0 +)
RATL	30	39	27	62	78	74 (+ +)	41 (0 +)	78 (+ +)	87 (+ +)
	60	48	23	52	55	72 (+ +)	36 (0 +)	69 (+ +)	71 (+ +)
	90	35	17	47	53	66 (+ +)	32 (0 +)	67 (+ +)	69 (+ +)
	120	34	19	50	62	65 (+ +)	34 (0 +)	68 (+ +)	73 (+ +)
RAVO	30	54	38	70	76	78 (+ +)	40 (0 +)	76 (+ +)	71 (+ +)
	60	58	37	64	74	76 (+ +)	36 (0 +)	72 (+ +)	67 (+ +)
	90	56	46	61	63	77 (+ +)	38 (0 +)	66 (+ +)	63 (+ +)
	120	52	35	50	64	69 (+ +)	34 (0 +)	64 (+ +)	60 (+ +)
RUDB	30	92	88	75	78	93 (+ +)	91 (+ +)	72 (+ +)	72 (+ +)
	60	87	90	48	67	91 (+ +)	92 (+ +)	60 (+ +)	64 (+ +)
	90	81	94	46	56	92 (+ +)	94 (+ +)	56 (+ +)	57 (+ +)
	120	88	92	44	56	89 (+ +)	92 (+ +)	55 (+ +)	56 (+ +)
RUDK	30	78	85	69	73	80 (+ +)	79 (+ +)	75 (+ +)	70 (+ +)
	60	76	80	54	68	75 (+ +)	74 (+ +)	64 (+ +)	64 (+ +)
	90	80	82	50	60	76 (+ +)	74 (+ +)	60 (+ +)	57 (+ +)
	120	72	80	52	50	74 (+ +)	73 (+ +)	56 (+ +)	53 (+ +)
RUDL	30	75	63	47	52	80 (+ +)	76 (+ +)	72 (+ +)	71 (+ +)
	60	73	64	32	47	76 (+ +)	73 (+ +)	64 (+ +)	63 (+ +)
	90	68	57	40	45	73 (+ +)	69 (+ +)	62 (+ +)	61 (+ +)
	120	71	66	41	42	74 (+ +)	70 (+ +)	59 (+ +)	57 (+ +)
RUDN	30	78	85	87	70	89 (0 0)	80 (0 0)	87 (+ +)	89 (+ +)
	60	52	86	57	71	38 (0 +)	78 (0 0)	69 (+ +)	71 (+ +)
	90	38	83	54	58	30 (0 +)	79 (0 0)	65 (+ +)	63 (+ +)
	120	31	84	57	54	27 (0 +)	80 (0 0