DOI: https://doi.org/10.37162/2618-9631-2022-2-6-52

УДК 551.509.313+551.509.324.2+551.508.85

 

 

Опыт пространственной верификации
радиолокационного наукастинга осадков:
определение и статистика объектов,
ситуаций и условных выборок

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

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

muravev@mecom.ru, boundel@mecom.ru

 

Проведен статистический анализ пространственных объектов и специальных ситуаций, определяемых для объектно-ориентированной верификации прогнозов областей осадков значительной и максимальной площади с привлечением методов теории экстремальных величин. Оценена степень влияния пропусков в точках полей и разных способов формирования пар прогноз-наблюдение на объемы и на статистические свойства выборок для пространственной верификации. С помощью пространственных квантильных функций и географических представлений продемонстрированы особенности композитных полей, построенных на основе данных около трех десятков радиолокаторов на Европейской территории России.

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

 

 

Expertise in spatial verification
of radar precipitation nowcasting:
identification and statistics of objects, situations
and conditional samples

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

Hydrometeorological Research Center of Russian Federation, Moscow, Russia

muravev@mecom.ru, boundel@mecom.ru

 

Statistical analysis was performed using methods of the extreme value theory for spatial objects and specified situations identified for object-oriented verification of precipitation regions with substantial and maximal areas. We made an estimation of the effect of missing values at field points and of different observation-forecast pairs construction on volumes and on statistical characteristics of samples retrieved for spatial verification purposes. We used spatial quantile functions and geographical representations in regular coordinates to illustrate particular aspects of composite fields built on about three dozen radars' data over the European territory of Russia.

Keywords: spatial forecast verification, radar precipitation nowcasting, extreme value theory, missing data, conditional verification sampling, spatial quantiles

 

Введение

На протяжении 2017–2020 гг. в теплый и холодный периоды года в Гидрометцентре России проводились испытания двух версий системы радиолокационного наукастинга осадков, которые обозначим по аналогии прогностическим моделям через STEPS-2014 и pySTEPS-2019. Особенности указанных моделей, условия и результаты испытаний изложены в нескольких работах авторов [13, 14, 20–23]. К теплому и холодному периодам отнесены по пять последовательных месяцев из интервалов май–сентябрь и ноябрь–март соответственно.

Первая версия модели STEPS-2014 испытывалась по территории Центрального федерального округа (ЦФО) на основе данных наблюдений девяти радиолокаторов ДМРЛ-С (р/л). Вторая версия pySTEPS-2019, усовершенствованная, испытывалась по Европейской территории России на основе данных 28 р/л (рис. 1).

 

а)                                                                                б)

Рис. 1. Область усовершенствованной системы наукастинга радиолокационных осадков: оперативная карта ЦАО и контуры характерных областей; сиреневый – начальная область анимации прогнозов (сайт ГМЦ); синий – объединенное поле круговых обзоров р/л (матрица 2151х1951); коричневый: условная карта для верификации (45°–65° с. ш., 30°–50° з. д.) (а); расположение ДМРЛ-С в условных координатах (28 р/л). Прямыми зелеными линиями на панели выделены ориентировочные направления на страны света (б).

Fig. 1. Area of improved radar precipitation nowcasting: operational map of the Central Federal district and the contours of characteristic areas; orange: initial area as in forecast animation (HMC site); blue: unified field of the circle radar area (matrix 2151x1951); brown: verification map in regular coordinates (45°-65°N, 30°-50°W) (a); Radar locations in special regular coordinates (28 radars). Straight green lines indicate the cardinal points (б).

 

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

Базовые, традиционные оценки качества наукастинга строились на временных рядах пар наблюдение-прогноз в узлах полей с последующим обобщением при помощи квантилей пространственного распределения этих оценок. Наряду с традиционными оценками рассчитывались показатели объектно-ориентированных и окрестных методов верификации [20, 31, 35, 36]. В рамках объектно-ориентированной верификации выделялись области осадков большой и максимальной площади, полная совокупность которых описывалась моделью "пиков над порогом" (POT, Peaks-over-Threshold) и анализировалась с помощью обобщенного распределения Парето из теории экстремальных величин [6–8, 17, 28, 46, 56]. Методы окрестной верификации были представлены показателем FSS (Fractions Skill Score) [48], который строится на пространственных долях областей наблюдаемых и прогнозируемых зон осадков.

Усовершенствованная система наукастинга была испытана и внедрена в оперативную практику в течение 2019–2020 гг. [23]. Здесь основное внимание уделено учету пропусков в данных наблюдений и разным способам формирования выборок пар наблюдение-прогноз. Информация о пропусках важна для оценок значимости характеристик качества, а испытание разных типов выборок позволяет учесть влияние на характеристики качества большого количества нулевых осадков и некоторого нижнего порога наличия осадков, учитывающего точность радиолокатора.

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

Указанные группы входят в классификацию методов пространственной верификации [38], в которую наряду с перечисленными тремя (объектно-ориентированные, окрестные и масштабно-разделительные) входит группа деформационных методов и иногда добавляется группа мер дальности [34]. Данная классификация не является строгой и здесь используется исключительно для отсылки на мировой опыт [3, 38, 39]. 

Уточним понятия областей, объектов и ситуаций. В полях осадков области выделяются, как обычно, пороговыми значениями. Ограниченные изолинией 1 мм/ч области будем называть объектами в соответствии с терминологией объектно-ориентированной верификации [31]. Площадь, или размер объекта оцениваются по количеству попадающих в него ячеек (или узлов), при таком определении связная одномерная «линия» всегда имеет ненулевую «площадь». Объекты значительной площади (значительные, или экстремальные объекты), составляющие выборку для приложения теории экстремумов, извлекаются из временной последовательности полей, которая отражает некоторую непрерывно развивающуюся синоптическую мезомасштабную ситуацию, или дождевую систему [1, 4, 32, 42]. Под точной непрерывной последовательностью полей понимается чередование их временных сроков в исходной 10-минутной дискретности; будем называть просто непрерывной такую последовательность, в которой, ввиду возможных аварий или нарушений информационного потока, допускаются один-два пропуска между соседними полями с возможностью интерполяции недостающих значений.

Пороговые значения определяются на основе документированного опыта и по результатам авторских экспериментов с полями осадков, а также с пакетами статистической обработки. Например, интенсивность 1 мм/ч для 10-минутных радиолокационных сумм является в климатологическом смысле достаточно высокой в средних широтах [1, 2, 11, 12], и в наших экспериментах такой порог обеспечивает не только заметные по площади объекты, но и временную непрерывность содержащих их полей в вышеуказанном смысле.

Объекты, ситуации и объемы выборок в точках полей осадков имеют важные статистические и географические свойства [31]; в частности, фундаментальные соотношения между пространственными и временными масштабами в метеорологических полях лежат в основе использованной модели наукастинга [13, 29, 51–53]. Заметим, что методы анализа пространственных характеристик областей осадков, включая экстремальные, давно применяются в гидрологии и содержатся в соответствующих нормативных документах ВМО [5, 57]. Описанию и анализу этих особенностей полей наблюдений и прогнозов в серии проведенных испытаний посвящена настоящая, первая статья предполагаемой серии. Качество прогнозов областей осадков значительной площади на основе обобщенного распределения Парето и качество прогнозов на основе отношения пространственных долей с помощью показателя FSS будут проанализированы во второй и третьей статьях предполагаемой серии соответственно; в этих работах областями испытаний являются круговые обзоры отдельных радиолокаторов. Результаты верификации прогнозов по полю, объединяющему обзоры 28 локаторов на Европейской территории России, планируется привести в четвертой статье.

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

 

1. Построение объединенного поля радиолокационных данных на Европейской территории России

Опишем область ответственности усовершенствованной схемы наукастинга 2020 г. и методы построения объединенного (или композитного) поля. На оперативной карте Центральной аэрологической обсерватории (ЦАО) (рис. 1а) эта область выделена синим полигоном, в нерегулярных координатах которого размещена сетка с 1951 точками по юго-восточной и северо-западной сторонам и 2151 точками по юго-западной и северо-восточной сторонам. Внутри этого полигона в настоящее время эксплуатируются двадцать восемь ДМРЛ-С: RUDO (Оренбург), RUDP (Петрозаводск), RUTD (Тамбов), RUDX (Архангельск), RAVO (Воейково), RAKT (Котлас), RUDB (Брянск), RUDV (Вологда), RUDG (Волгоград), RATL (Тула), RAKD (Краснодар), RUDT (Ставрополь), RUWJ (Валдай), RAMI (Миллерово), RAVN (Внуково), RAKU (Курск), RUDI (Ижевск), RAYL (Элиста), RUDL (Смоленск), RASM (Самара), RABG (Белгород), RUDM (Мин. Воды), RUDK (Кострома), RUDZ (Казань), RAVL (Вел. Луки), RUDU (Уфа), RUDN (Ниж. Новгород), RAKW (Киров). В предыдущей схеме наукастинга и в испытаниях 2017–2019 гг. были задействованы девять р/л, покрывающих в основном территорию ЦФО: RAVO, RUDB, RATL, RUWJ, RAVN, RAKU, RUDL, RUDK, RUDN.

Объединенные поля наблюдений и прогнозов, как и поля показателей качества, отображаются на условную географическую карту (рис. 1б), на которой сохранена структура расположения локаторов, а географическая сеть задана в регулярных широтно-долготных координатах 45–65° и 30–50° с подходящими инкрементами для вложения матрицы размером 2151×1951.

Применение условной картографии оправдано тем, что в рабочей версии системы Grads существует только два способа географического представления данных – в виде регулярной географической сетки или с помощью точечных координат (в нашем случае в количестве порядка четырех миллионов). Очевидно, в последнем случае и составление координат в текстовом формате, и выведение в графику с помощью внутренних процедур сглаживания превращаются в трудно исполнимое предприятие. К тому же детальная географическая картина на данном этапе анализа не требуется. Названия пунктов на рис. 1б в центрах кругов соответствуют спецификации оперативной карты радиолокационной сети ЦАО, прямыми линиями выделены ориентировочные направления на страны света.

Опишем процедуру формирования композитного поля. Входные композиты размером 2151×1951, заполненные вначале признаками отсутствия данных (числом -999), замещаются в соответствующих координатах оценками интенсивности осадков, которые формируются в ЦАО в виде квадратных матриц размером 504×504 с вложенными круговыми зонами обзора. В областях перекрытия зон обзора локаторов значение интенсивности рассчитывается с помощью линейной весовой формулы с приоритетом ближайшего локатора.

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

Выходные поля модели не содержат признаков отсутствия данных, кроме тех случаев, когда движение начинается от края видимой зоны. Точки с константами отсутствия в исходных радиолокационных полях переносятся на каждом шаге прогноза в выходное поле модели. Несомненно, полезное свойство усовершенствованной схемы состоит в том, что адвекция точек поля с числовым значением позволяет заполнять осадками некоторые точки, которые исходно были "пустыми" [54].

 

2. Выборки объектов значительной площади: ситуации с пиком

В нашем подходе рассматривается оценка способности модели воспроизводить статистические характеристики тех областей осадков, площадь которых либо максимальна, либо превышает некоторый пороговый уровень. Как известно, распределения экстремальных величин описываются теоремой Фишера – Типпета – Гнеденко (первая теорема теории экстремальных величин, ТЭВ), а распределения превышений порогового уровня описываются теоремой Пикандса – Балкема – де Хаана (вторая теорема ТЭВ). В рамках каждой из этих моделей существует по три основных закона притяжения, т. е. три класса распределений, к которым асимптотически стремятся распределения некоторым образом нормированных экстремумов (см., например, монографии [9, 10, 28, 33, 46, 47]). Так как исходные данные являются временными последовательностями полей, то статистический анализ может рассматриваться в рамках теории выбросов ("максимумов", "пиков", "превышений") случайных процессов (см., например, монографии [6, 7, 15, 17, 18, 26]). 

Проблемы моделирования на основе этих двух теорем ТЭВ относятся к "классической теории экстремальных величин", считающейся в некотором смысле завершенной [7]. Однако статистический анализ экстремальных величин по конкретному набору данных (статистический вывод) оказывается намного более сложным и далеким от завершения, что отмечается во всех литературных источниках по ТЭВ, затрагивающих проблемы приложений теорем к практике. Эти сложности относятся в первую очередь к проверке условий используемых теорем ТЭВ. Доказать корректность подготовленной выборки (что равносильно проверке соответствия выборки условиям соответствующей теоремы) не только затруднительно, но скорее всего невозможно: "вопрос о проверке выполнения в реальности условий теоремы чаще всего не может быть не только решен, но даже и поставлен" [27]. При моделировании пиков над порогом вышеуказанная проверка превращается в корректный выбор порога, а так как для анализа экстремальных значений временных рядов выборка формируется с учетом временной локальности пиков, то приходится дополнительно определять подходящие временные интервалы – носители искомых экстремумов. Для подтверждения условий второй теоремы ТЭВ мы учитываем физические и эмпирические представления о характерных пространственно-временных особенностях мезомасштабных осадков [1, 4, 12, 31, 32, 42], а также пуассоновское распределение моментов времени появления экстремумов над выбранными порогами [8, 15, 18, 24, 25, 30, 41, 44, 46, 56].

Конкретизируем некоторые понятия и термины. Будем называть носители искомых пиков ситуациями с пиком, или ситуациями, со следующим определением – это последовательность полей осадков с интенсивностью не менее 1 мм/ч хотя бы в одной точке в каждом поле и с наличием в этой последовательности хотя бы одного поля с областью сплошных осадков, площадь которой превышает заданное пороговое значение. В данном случае сплошность, или связность области являются характеристиками пространственного осреднения. Порог 1 мм выбран намеренно таким же, как и порог для идентификации объектов.

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

Задача комплектации значительных объектов в выборки для верификации с помощью распределения Парето решалась экспериментально на основе анализа гистограмм распределения размеров всех объектов с особым вниманием к правым хвостам этих распределений. Процедуры выделения ситуаций применялись к последовательностям полей наукастинга всех заблаговременностей от 10 до 120 минут. Пороги превышения определялись также экспериментально по набору 625, 900, 1225 и 1600 точек поля и с предпочтением таких значений, которые приводят к наиболее устойчивым оценкам параметров распределения Парето.

 

3. Особенности определения и выделения объектов

Выделение объектов производится с помощью функции FeatureFinder библиотеки SpatialVx [37] из хранилища языка R на основе следующих параметров:  

– порог идентификации – 1 мм/ч; 

– параметр свертки (радиус осреднения) – 9 узлов сетки (~18 км); 

– минимальные размеры (они же пороговые значения размера, или параметры положения распределения Парето, или параметры Парето) – 625=25´25, 900, 1225, 1600=40´40 узлов сетки (линейный размер одной ячейки равен 2 км).

Подчеркнем еще раз, объект выделяется после пространственного сглаживания поля.

Рассмотрим пример функционирования модели STEPS и пакета SpatialVx на осадках локатора RAKU (Курск). Обратим внимание только на общую форму объектов без указания значений интенсивности. В качестве начальных данных модели STEPS берутся два последовательных поля отражаемости размером 504´504 км, при этом радиолокационные наблюдения заполняют вписанный круг обзора, вне которого размещаются константы отсутствия. Для использования в модели эти поля интерполировались в сетку 256×256 с разрешением около 2 км, а в процессе преобразования отражаемости в интенсивность осадков константы отсутствия вне круга обзора заменялись небольшими положительными числами порядка 10-3. Указанные манипуляции требовались по условиям функционирования модели STEPS, в которой двумерное преобразование Фурье применялось к полному числовому полю данных, а линейный размер поля должен был выражаться степенью двойки.   

Расчет характеристик объектов и графическое представление самих объектов с помощью пакета SpatialVx начинается с прогнозов на первый срок 10 минут (рис. 2а). В этом примере порог равен 1225 точек (эквивалент квадрата со стороной 70 км); хорошо видно, что радиус осреднения в 9 узлов и минимальный порог площади устранили множество мелких несущественных объектов. Построение пар сравниваемых объектов производится по минимальному расстоянию между центроидами полей наблюдения и прогноза (т. е. формально рассчитанными "центрами тяжести объектов"); спаренные объекты отмечаются одинаковым цветом.

Как видно по первой строке панелей (рис. 2а), уже на следующем шаге регрессионной схемы происходят несущественные, но заметные изменения выделенного объекта прогнозирования. К сроку 60 мин сам объект (внутри изолинии 1 мм) теряет в площади, но мало изменяет ориентацию, хотя объект в поле наблюдения развернулся примерно на 10–25о против часовой стрелки с усилением полосы интенсивности в южной части, направленной на CВ. К концу прогностического интервала (рис. 2в) объект наблюдения разбивается на два, а прогностический объект, заметно отставая во вращении, переключается алгоритмом спаривания на "отделившийся аппендикс" – верификация прогноза первоначального объекта теряет смысл. К этому сроку прогноза из атрибутов объектов, по которым еще можно оценивать качество прогноза, осталась, по-видимому, только максимальная площадь объекта в полном поле, а поскольку подобная "потеря пары" встречается довольно часто, то выбор максимальной площади объекта в качестве величины для верификации получает дополнительный кредит.

 

а)

б)

в)

Рис. 2. Эволюция полей и объектов в наблюдении (RAKU, 1-й и 3-й столбцы) и в прогнозе (STEPS, 2-й и 4-й столбцы) через 10 мин (а), 60 мин (б) и 150 мин (в). Затенённые столбцы слева – эволюция полей осадков в исходном разрешении (сторона ячейки около 2 км), светлые столбцы справа – эволюция выделенных объектов с заданными параметрами порога интенсивности осадков (1 мм/ч), радиусом осреднения 9 точек поля (18 км) и минимальным размером 35×35. Красным цветом в двух правых столбцах выделены спаренные объекты, для которых далее рассчитываются характеристики ошибок. Синий объект в наблюдении (в) – объект, оставшийся к концу прогностического интервала "без пары".  

Fig. 2. Evolution of fields and objects in the observations (RAKU radar, the 1st and 3rd columns) and in the forecast (STEPS, 2nd and 4th columns) in 10 min (a), 60 min (b), and 150 min (c). Shadowed columns on the left represent the evolution of precipitation fields in the original resolution (grid cell side of about 2 km), the light-colored columns on the right are the evolution of the objects determined with the following parameters: precipitation threshold (1 mm/h), convolution radius of 9 grid points (18 km), and the minimum size of 35×35 (70×70 km square). Red in the columns on the right indicates matched objects, for which the scores are calculated. The blue object in (c) is "unpaired".

 

 

3.1. Статистические характеристики ситуаций и объектов

Будем предполагать, что в последовательности полей, организованных в 10-минутной дискретности и обладающих объектами достаточно большого размера, можно возместить отсутствие одного поля или двух соседних полей интерполированием, пользуясь инерцией физического процесса: можно сказать, ситуация сшивается по интервалу сцепления в 20 мин. Максимальный размер объекта в каждой ситуации с пиком послужит выборочным значением для оценки распределения Парето с параметром положения из указанного выше набора 625, 900, 1225 и 1600 точек.

Рассмотрим временной ход площади максимального объекта в последовательных полях обзора локатора Курск в период 5–19 мая 2017 г. (рис. 3); на график наносятся только значения не меньше 625 точек. Сроки с пиками отмечены синими линиями; предполагается, что они случайно разбросаны по оси времени в соответствии с простым потоком Пуассона. Отметим несколько особенностей на рис. 3. Как оказалось, в интервале 10.50–12.20 ч ВСВ наблюдения полностью отсутствовали, поэтому пик 0508_1020 следовало бы аннулировать, что, однако, нельзя было исполнить в режиме автоматизированной обработки. Наконец, серия значительных площадей в конце выделенного интервала (до срока 0519_1930) относится к ситуации, выходящей за пределы интервала данного рисунка. 

Рис. 3. Площади объектов в последовательных полях наблюдений р/л Курск в период с 5 по 19 мая 2017 г., превышающие порог 625 точек (красные точки). Синие линии соответствуют срокам и максимальным площадям среди объектов, принадлежащих отдельной "ситуации". На оси времени представлены все сроки наблюдений, интервалы с отсутствием точек указывают либо на пропуски, либо на отсутствие объектов, имеющих размер не менее 625 точек.

Fig. 3. Areas of objects not less than the threshold of 625 grid points (red points) in consecutive observation fields of the Kursk radar from 09 May to 19 May 2017. Blue lines correspond to times and the maximum areas among the objects belonging to a situation. The time axis shows all the observation times, intervals without points indicate either missing data or absence of the objects not less than 625 grid points.

 

 

Точная информация об объектах в нескольких последовательных сроках выделенного интервала показана в табл. 1. В частности, видно наличие полей с набором из нескольких подходящих объектов, включая объекты размером в одну точку.

 

Таблица 1. Пример ситуации "0507_0910 – 0507_1810" с указанием в полях наблюдения локатора Курск площадей всех выделенных объектов, максимальной площади среди объектов одного поля и "пиковой" площади всей ситуации, используемой в выборке пиков над порогом. Отобраны максимальные площади размером не менее 625 точек

Table 1. Example of the situation "0507_0910 – 0507_1810" for the Kursk radar, with the areas of all identified objects, maximum area among the objects of one field, and the “peak” area of the situation used in the peaks-over-threshold sample. Maximum areas of not less than 625 grid points are used

Срок

Набор площадей
всех объектов поля

Суммарная
площадь

Максим.
площадь

ПИК
"ситуации"

0507_0830

0

 

 

0507_0840

0

 

 

0507_0850

0

 

 

0507_0900

0

 

 

0507_0910         215

215

 

 

0507_0920        

0

 

 

0507_0930         257

257

 

 

0507_0940         301

301

 

 

0507_0950         308

308

 

 

0507_1000         422

422

 

 

0507_1010         912

912

912

 

0507_1020         769  32  27  89

917

769

 

0507_1030         892  127  18  27  72

1136

892

 

0507_1040         1086  188  1

1275

1086

 

0507_1050         1038  184

1222

1038

 

0507_1100         903  235  24

1162

903

 

0507_1110         1237  44  1  16

1298

1237

 

0507_1120         1350  19  29  66

1464

1350

 

0507_1130         1376  568  10  31  18

2003

1376

 

0507_1140         1377  749  21  17  8

2172

1377

 

0507_1150         1666  18  12  10

1706

1666

1666

0507_1200         1637  6

1643

1637

 

0507_1210         452  929

1381

929

 

0507_1220         421  188  55  20  653

1337

653

 

...

...

 

 

0507_1710         491

491

 

 

0507_1720         470

470

 

 

0507_1730         505

505

 

 

0507_1740         556

556

 

 

0507_1750         569

569

 

 

0507_1800         105

105

 

 

0507_1810          99

99

 

 

0507_1820

0

 

 

0507_1830

0

 

 

Соберем в один рисунок несколько характерных полей наблюдений, содержащих пики, и соответствующие прогнозы на эти сроки со стартом на час раньше (рис. 4).

 

 а)

 б)

 в)

 г)

 д)

Рис. 4. Левый столбец панелей: поля радиолокационных осадков, полученных по данным наблюдений локатора Курск и содержащих объекты пиковой площади, выделенные на рис. 3. синими линиями. Площади объектов (в количестве точек): 8516 (а); 3747 (б); 7895 (в); 5337 (г); 8098 (д). Правый столбец панелей: прогноз на эти сроки наблюдений, выпущенный часом раньше.

Fig. 4. Left column: the Kursk radar precipitation fields containing the objects with peak areas indicated with the blue lines in fig. 3. Areas of objects in grid points: 8516 (а); 3747 (б); 7895 (в); 5337 (г); 8098 (д). Right column: forecast for this time issued one hour earlier.

 

 

Необходимо иметь в виду, что поле наблюдения состоит из фактических значений интенсивности осадков только внутри вписанного круга, границы которого заметно проявляются при наличии осадков (например, левые панели а, б, в, г); в то же время "пустоты" в полях прогноза могут заполняться прогностическими данными (правые панели а, г). Добавим, что на рис. 3 фактические поля наблюдений представлены в разрешении 1 км, а поля соответствующих прогнозов – в разрешении 2 км.

В табл. 2 и 3 собраны суммарные характеристики всех выделенных ситуаций с учетом количества и размеров объектов в полях обзора локатора RAKU (Курск) и в соответствующих прогностических полях за пять месяцев теплого периода 2017 г. Табл. 2 и 3 отличаются лишь минимальным размером идентифицируемых объектов – не менее 625 и не менее 1600 точек поля, соответственно. Пики ситуаций – это значения площади объектов в столбце max.    

Сделаем несколько замечаний. Первое замечание относится к методу выделения ситуаций при разных порогах. Сравним первые строчки первых частей обеих таблиц, относящихся к радарным наблюдениям: ситуация длительностью 550 минут (9 часов 10 минут) содержит 19 объектов размером не менее 625 точек и всего 2 объекта размером не менее 1600 точек. Так как несколько объектов могут «населять» одно поле, то их формально просуммированное по каждому полю количество (в столбце objs) может существенно превышать количество полей в ситуации. Однако количество синоптических объектов мезомасштабных размеров, изменяющихся в пространстве и во времени на протяжении одной ситуации, как правило, на порядки меньше общего количества ситуаций (см. по этому поводу обсуждение эволюции объектов во времени в статье [32]). А ввиду вложенности ситуаций друг в друга "по образцу матрёшки" (см. табл. 1), объекты размером не меньше 1600 точек находятся в выборке объектов размером не менее 625 точек – в набор пиковых значений в обоих случаях попадает одно и то же значение 1666. Дополнительно отметим, что во втором случае – при двух выделенных объектах – все квартили оказываются различными из-за линейной интерполяции, автоматически проводимой на не менее двух значениях.

Второе замечание касается количества выделенных ситуаций в прогнозах в сравнении с наблюдениями. Автоматическое выделение ситуаций нуждается иногда в корректировке из-за сбоев в технологии прогнозирования. Так, в период первого испытания технологии наукастинга в 2017 г. потери в прогнозах составляли до 10–12 % от максимально возможного количества в 22032 10-минутных интервалов [21, табл. 1]. Такие потери означают пробелы в непрерывной последовательности полей, приводящие иногда к короткоживущим ситуациям длительности от 10 до 30 мин. С учетом таких потерей можно сделать вывод о примерно 20%-м росте количества выделенных ситуаций в прогностических полях при прогнозе на 60 мин для обоих порогов 625 и 1600 точек.


 

 

Таблица 2. Характеристики выделенных «ситуаций» в ряду наблюдений и наукастинга по данным радиолокатора Курск в теплый период года (май-сентябрь 2017 г.). Размер идентифицируемого объекта – не менее 625 точек поля

Table 2. Characteristics of indentified situations in observations and nowcasting data series from the Kursk radar during the warm period (May-September 2017). Оbject size not less than 625 grid points

sit

min

q25

med

q75

max

objs

drt

d_time_start

d_time_stop

Радиолокатор КУРСК

1

653

754

912

1294

1666

19

550

2017-05-07 09:10

2017-05-07 18:10

2

666

714

907

924

1022

5

150

2017-05-08 00:50

2017-05-08 03:10

3

893

1529

1828

2770

3292

11

140

2017-05-08 08:30

2017-05-08 10:40

4

678

2112

2895

5790

8516

67

750

2017-05-08 12:30

2017-05-09 00:50

5

647

1412

1720

2460

3747

56

720

2017-05-09 20:50

2017-05-10 08:40

6

645

1694

2910

3886

7087

141

1470

2017-05-11 14:00

2017-05-12 14:20

7

634

1202

2444

4677

7895

111

1400

2017-05-12 15:10

2017-05-13 14:20

8

628

1184

1676

2110

2600

49

720

2017-05-13 16:10

2017-05-14 04:00

9

629

1210

1753

3032

5337

108

1130

2017-05-16 08:20

2017-05-17 03:00

10

628

954

1584

2328

8098

386

3970

2017-05-17 04:40

2017-05-19 22:40

. . . . . . . . . . .

83

638

768

899

985

1071

3

300

2017-09-18 10:20

2017-09-18 15:10

84

629

676

704

747

1040

10

270

2017-09-18 17:20

2017-09-18 21:40

85

684

1396

2974

5066

7720

140

1670

2017-09-20 00:00

2017-09-21 03:40

86

627

855

1090

1188

1339

15

420

2017-09-24 05:40

2017-09-24 12:30

Прогностическая система STEPS

1

680

1146

1426

1757

2296

30

760

2017-05-07 09:20

2017-05-07 21:50

2

685

753

821

900

978

3

300

2017-05-07 22:30

2017-05-08 03:20

3

634

1142

1914

2284

3187

8

130

2017-05-08 08:40

2017-05-08 10:40

4

655

2056

4079

6048

9031

67

750

2017-05-08 13:00

2017-05-09 01:20

5

646

662

672

705

794

4

340

2017-05-09 14:50

2017-05-09 20:20

6

1113

1732

2468

3153

4380

14

180

2017-05-09 21:10

2017-05-10 00:00

7

684

1536

1964

2556

3580

42

620

2017-05-10 00:30

2017-05-10 10:40

8

708

722

767

911

918

5

200

2017-05-10 15:00

2017-05-10 18:10

9

628

1418

2560

3814

7514

140

1480

2017-05-11 14:20

2017-05-12 14:50

10

659

2407

3780

6488

8523

45

490

2017-05-12 15:50

2017-05-12 23:50

. . . . . . . . . . .

115

628

1341

2537

4740

6824

119

1550

2017-09-19 22:10 

2017-09-20 23:50

116

654

1399

1635

2763

3667

17

290

2017-09-21 00:20 

2017-09-21 05:00

117

637

638

640

664

688

3

210

2017-09-23 03:20 

2017-09-23 06:40

118

756

913

977

1131

1598

14

420

2017-09-24 06:00 

2017-09-24 12:50

Примечание. Интервал сцепления равен 20 мин. Минимальный размер связного объекта составляет 625 точек поля. Экстремумы и оценки квартилей размеров всех объектов ситуации (в количестве внутренних точек) размещены в столбах min, q25 (25%-й квартиль), med (медиана), q75(75%-й квартиль) и max. В столбцах objs и drt находятся количество пиков (объектов размера не меньше 625 точек в данной ситуации) и длительность ситуации (в минутах) соответственно. Два последних столбца содержат начальную и конечную дату и время ситуации (включительно).  

 

 

Таблица 3. Характеристики выделенных «ситуаций» в ряду наблюдений и наукастинга по данным радиолокатора Курск в теплый период года (май-сентябрь 2017 г.). Размер идентифицируемого объектане менее 1600 точек поля

Table 3. Characteristics of indentified situations in observations and nowcasting data series from the Kursk radar during the warm period (May-September 2017). Оbject size not less than 1600 grid points

sit

min

q25

med

q75

max

objs

drt

d_time_start

d_time_stop

Радиолокатор КУРСК

1

1637

1644

1652

1659

1666

2

550

2017-05-07 09:10

2017-05-07 18:10

2

1623

1825

2432

2997

3292

8

140

2017-05-08 08:30

2017-05-08 10:40

3

1635

2176

3026

6257

8516

62

750

2017-05-08 12:30

2017-05-09 00:50

4

1651

1886

2415

2656

3747

33

720

2017-05-09 20:50

2017-05-10 08:40

5

1609

2182

2977

4256

7087

111

1470

2017-05-11 14:00

2017-05-12 14:20

6

1644

2316

3336

6098

7895

78

1400

2017-05-12 15:10

2017-05-13 14:20

7

1608

1789

2088

2256

2600

28

720

2017-05-13 16:10

2017-05-14 04:00

8

1600

1805

2325

3535

5337

70

1130

2017-05-16 08:20

2017-05-17 03:00

9

1605

1863

2362

3580

8098

190

3970

2017-05-17 04:40

2017-05-19 22:40

10

1607

1749

2170

2697

3815

71

790

2017-05-21 09:30

2017-05-21 22:30

. . . . . . . . . . .

42

1621

1801

2056

2324

4489

28

830

2017-09-05 07:50

2017-09-05 21:30

43

1602

1638

1658

1770

1844

8

590

2017-09-10 00:50

2017-09-10 10:30

44

1603

1603

1603

1603

1603

1

300

2017-09-13 11:10

2017-09-13 16:00

45

1602

2648

3572

5756

7720

103

1670

2017-09-20 00:00

2017-09-21 03:40

Прогностическая система STEPS

1

1643

1713

1904

1981

2296

12

760

2017-05-07 09:20

2017-05-07 21:50

2

1793

2034

2249

2391

3187

5

130

2017-05-08 08:40

2017-05-08 10:40

3

1601

2396

4501

6270

9031

59

750

2017-05-08 13:00

2017-05-09 01:20

4

1673

2188

2947

3444

4380

11

180

2017-05-09 21:10

2017-05-10 00:00

5

1616

2019

2291

2860

3580

28

620

2017-05-10 00:30

2017-05-10 10:40

6

1629

2448

3222

4135

7514

98

1480

2017-05-11 14:20

2017-05-12 14:50

7

1823

3324

5638

6702

8523

35

490

2017-05-12 15:50

2017-05-12 23:50

8

1620

1952

2492

4519

6871

42

920

2017-05-13 00:20

2017-05-13 15:30

9

1604

1906

2165

2482

3264

26

950

2017-05-13 16:40

2017-05-14 08:20

10

1631

2044

2704

4121

5841

54

930

2017-05-16 08:30

2017-05-16 23:50

. . . . . . . . . . .

61

1604

1849

2124

2463

3364

30

730

2017-09-05 08:20

2017-09-05 20:20

62

1642

1719

1838

1970

2215

10

600

2017-09-10 01:50

2017-09-10 11:40

63

1605

2684

4066

5366

6824

78

1550

2017-09-19 22:10

2017-09-20 23:50

64

1635

2058

2763

2949

3667

9

290

2017-09-21 00:20

2017-09-21 05:00

Примечание. Интервал сцепления равен 20 мин. Минимальный размер связного объекта составляет 1600 точек поля. Экстремумы и оценки квартилей размеров всех объектов ситуации (в количестве внутренних точек) размещены в столбах min, q25 (25%-й квартиль), med (медиана), q75(75%-й квартиль) и max. В столбцах objs и drt находятся количество пиков (объектов размера не меньше 1600 точек в данной ситуации) и длительность ситуации (в минутах) соответственно. Два последних столбца содержат начальную и конечную дату и время ситуации (включительно).  


 

 

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

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

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

 

Таблица 4. Диапазоны изменений минимальных (min) и максимальных (max) размеров объектов в выделенных ситуациях. Минимальный размер связного объекта – 625 точек поля

Table 4. Range of minimal (min) and maximal (max) object sizes in the identified situations, Minimal object size is 625 grid points

Поля

Колич.
ситуаций

Характеристики выделенных ситуаций и диапазоны
 их изменений

min

max

objs

drt

мин

макс

мин

макс

мин

макс

мин

макс

RAKU

86

625

893

625

17989

1

386

90

3970

STEPS-30

117

625

1242

637

18282

1

183

10

1930

STEPS-60

118

625

1113

629

17601

1

185

10

2350

Примечание. objs – количество объектов в полях; drt – продолжительность ситуации в минутах. Интервал сцепления равен 20 мин.

 

 

 

В табл. 5 и табл. 6 представлено количество выделенных ситуаций с параметром сцепления 20 мин по областям обзора всех локаторов (кроме р/л Внуково) для заблаговременностей 30, 60, 90 и 120 минут и для двух пороговых площадей объектов – 625 и 1600 точек.

 

Таблица 5. Количество выделенных ситуаций в последовательностях полей наблюдений радиолокаторов и прогнозов разной заблаговременности с помощью системы STEPS в теплый период 2017 г. (май-сентябрь). Минимальный размер выделяемого связного «объекта» в поле равен 625 точкам. Параметр сцепления – 20 минут

Table 5. Number of the situations in consecutive fields of radar observations and STEPS nowcasts for different lead times in the warm period (May-September). The minimal object size is 625 grid points. Tie interval is 20 min – two consecutive fields missing in a sequence are allowed

Идентификатор
р/л

Наблюдения
р/л

Прогностические поля:
заблаговременность (мин)

30

60

90

120

RAKU

86

117

118

121

123

RATL

79

117

117

128

132

RAVO

89

133

141

145

157

RUDB

76

90

96

99

99

RUDK

96

130

142

145

155

RUDL

92

119

123

131

130

RUDN

84

122

129

141

147

RUWJ

85

119

128

138

152

 

 

Таблица 6. То же, что и в табл. 5, кроме минимального размера выделяемого связного «объекта» в поле: здесь он равен 1600 точкам 

Table 6. As in Table 5, but for the minimal object size of 1600 grid points

Идентификатор
р/л

Наблюдения
р/л

Прогностические поля:
заблаговременность (мин)

30

60

90

120

RAKU

45

65

64

64

67

RATL

56

71

73

76

76

RAVO

64

81

83

85

95

RUDB

51

57

57

59

60

RUDK

65

79

88

94

96

RUDL

60

79

83

85

92

RUDN

60

87

89

93

95

RUWJ

85

119

128

138

152

 

 

Сопоставим основные характеристики ситуаций и объектов с гистограммами распределения всех размеров объектов, выделенных порогом 1 мм/ч. Ограничимся гистограммами для данных локатора Курск (рис. 5), ввиду того, что для остальных локаторов картина сходная.

 

 

Рис. 5. Гистограммы распределения размеров всех объектов, выделенных порогом 1 мм/ч  в полях интенсивности осадков  радиолокатора Курск. Верхний ряд панелей – данные по теплому периоду, нижний – по холодному. Красный цвет – количество, синий – частоты. В названии панелей указаны количество связных объектов (cases), количество градаций гистограммы (breaks +1, в левом столбце) и основные характеристики распределения размеров (квартили и экстремальные значения). Для удобства графического представления размер объекта на оси абсцисс поделен на 1000.

Fig. 5. Histograms of the distribution of sizes of all objects identified using the threshold of 1 mm/h in the Kursk radar fields. The upper row is for the warm period, the lower row is for the cold period. Red indicates the number, blue, the frequencies. The number of contiguous objects (cases), histograms categories (breaks + 1, left column), and the basic distribution characteristics (quartiles and extremes) are given above the panels. The object size in abscissa axis is divided by 1000 for convenience.

 

 

Количество всех выделенных порогом 1 мм/ч объектов в теплый период примерно вдвое превышает аналогичное количество в зимний период (15403 против 7384, рис. 5, левый столбец), однако характер плотности распределения в основном одинаков (рис. 5, правый столбец). Сопоставим общие статистические характеристики распределения размеров объектов по данным теплого и холодного периода: медиана 317 против 262, 75%-й квартиль 1133 против 1080 и максимальное значение 17989 против 14228. Переход к линейным масштабам объектов (извлечением квадратного корня) позволяет сопоставить полученные результаты с данными [31].  

Распределение указанных пиковых площадей моделируется в дальнейшем обобщенным распределением Парето с помощью пакета extRemes языка R [40].

 

3.2. Пуассоновость распределения пиков на временной шкале

Известно, что обобщенное распределение Парето возникает как предельное распределение "пиков над порогом", но не всегда указывается, что выборка пиков предполагается взятой из обобщенного распределения экстремальных величин. Более того, если пики расположены на временной оси, как на примере рис. 3, то для выполнения условий первой теоремы об экстремумах достаточно соблюдения пуассоновости их распределения. Это предполагает, что отсчеты времени с пиками распределены как точечный (или простейший) пуассоновский поток событий [15, 17, 18, 44, 55], определяемый тремя условиями:

1) стационарность – вероятность появления k событий в любом промежутке времени длины s зависит только от k и s и не зависит от начала отсчета данного промежутка;

2) отсутствие последействия – появление k событий в промежутке времени не зависит от того, появлялись ли эти события ранее;

3) ординарность – появление не менее двух событий за малый промежуток времени "практически" невозможно.

Всем этим требованиям удовлетворяет следующее аналитическое определение потока Пуассона. Обозначим через λ среднее число событий в единицу времени (интенсивность потока). Если вероятность появления k событий за время t описывается формулой

 ,

то условия 1 3 автоматически выполняются.

Пуассоновское распределение нетрудно смоделировать с помощью серий независимых псевдослучайных чисел, воспользовавшись обратной функцией показательного распределения для сумм случайных величин (например, [16, стр. 60]). С помощью пакета poisson (из архива R) соответствие потоку Пуассона указанных на рис. 3 сроков было проверено на ансамбле 100 реализаций с последующими оценками средних, медианных, наилучших и наихудших последовательностей точек на оси времени. Интенсивность потока рассчитывается по всему периоду верификации, состоящему округленно из 22000 случаев. Нумерацию сроков увеличим примерно на 200 (предыдущие сутки не имели пиков) и положим длину интервала рис. 3 округленно равной 2000. В архиве наблюдений локатора Курск в теплый период 2017 г. выделено 86 пиков (или ситуаций), т. е. интенсивность потока равна λ = 0.0039, а значение λt = 7.802 для интервала рис. 3 дает ожидаемые около восьми пиков. Рис. 6 содержит все реализации точечного процесса Пуассона с параметрами от 1 до 10 для целых значений распределения Пуассона. На веер построенного ансамбля нанесены штриховой линией квантили 2.5 % и 97.5 % и желтой линией – ожидаемое среднее значение.  

Рис. 6. Ансамбль 100 реализаций потока Пуассона для моментов времени, в которые наблюдались пиковые площади выделенных объектов. Синими линиями отмечены: на оси времени условная длина интервала рис. 3 (2000 сроков), на оси количества пиков – параметр распределения Пуассона, соответствующий ожидаемому количеству пиков на данном интервале. Желтая диагональ соответствует ожидаемому количеству пиков в зависимости от длины интервала. Красные штриховые линии выделяют квантили 2.5 % и 97.5 % номеров соответствующих отсчетов времени; например, верхняя линия квантилей располагает до 8 пиков в первой тысяче сроков, в то время как нижняя линия располагает во всем интервале рис. 3 не более 4 пиков.

Fig. 6. Ensemble of 100 Poisson flow realizations for the time moments of peak object areas. Blue lines  on X-axis indicate conditional interval duration of Fig. 3 (2000 times), the axis of peak numbers shows the Poisson distribution parameter corresponding to the expected peak number on this interval. The yellow diagonal corresponds to the expected peak number depending on the interval length. The red dotted lines indicate quantiles 2.5% and 97.5% of the numbers of corresponding time counts; e.g., the upper line of quantiles places up to 8 peaks in the first thousand of times, while the lower line places not more than 4 peaks in the whole interval of fig. 3.

 

 

Простейшие характеристики данного ансамбля расположены в табл. 7. Видно, что номера сроков пиковых площадей по рис. 6. входят во внутреннюю часть веера реализаций; наименьшее среднее абсолютное расстояние от номеров отсчетов 100 реализаций указывает на ближайшую реализацию под номером 40. Однако и наблюдаемые сроки, и сроки ближайшей модельной реализации занижены по отношению к медианной и средней реализациям, что, возможно, объясняется наличием мнимых и не устраненных пиков на рис. 3, а также случайным характеров самого процесса генерирования пуассоновского потока. Приведенный пример проверки пуссоновости процесса может быть легко обобщен на полные интервалы испытаний.

 

 

Таблица 7. Номера сроков наблюдений (xobs) и квантилей распределения ансамбля точек отсчета модельного потока Пуассона на интервале до 2000 начальных сроков. Интенсивность потока Пуассона рассчитана для 86 пиковых площадей, не меньших 625 точек. Статистические характеристики модельного ансамбля - min, med, mean, max; ближайшая к наблюдениям реализация сроков - best

Table 7. Numbers of the times (xobs) and of the quantiles of ensemble distribution of the starting points of the modeled Poisson flow in the interval of up to 2000 initial times. The Poisson flow intensity is computed for 86 peak areas not less than 625 grid points. The statistical characteristics of the model ensemble are min, med, mean, max; best is the member closest to the observations

 

Пуассоновская реализация десяти последовательных сроков
 наблюдений в интервале 1–2000

      1         2         3        4          5         6        7         8         9        10

xobs

    217     305     352     380     573     898  1002   1149   1565   1661

best

    247     281     379     491     708     781    932   1063   1330   1680

min

        5       16     116     202     315     340     485    802     918   1071

med

    180     412     656     871   1108   1292   1577   1942   2184   2411

mean

    244     476     744     969   1261   1493   1778   2047   2291   2538

max

  1161   1825   2419   2691   3328   3454   4024   4148   4597   4680

 

 

4. Выделение полей с объектами значительной площади для оценок FSS

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

Для идентификации объектов в обоих периодах испытаний оставлен один порог (1 мм/ч). Локатор Внуково был исключен из-за незначительного количества выделенных ситуаций в холодный период. Количество накопленных полей с объектами площади не менее 1600 точек указано в табл. 8. Так как выборки составлялись из полей, принадлежащих зачастую одной ситуации, то оказался нарушенным принцип статистической независимости элементов выборки. Этот принцип был отчасти соблюден случайными отборами полей в количестве, не превышающем 200, что имеет порядок количества выделенных ситуаций в теплый период года.  

 

Таблица 8. Общее количество прогностических полей (сроков) с наличием объектов, выделенных после пространственного сглаживания с радиусом 9 узлов с помощью изолинии 1 мм/ч и имеющих размер не менее 1600 точек. Теплый период май-сентябрь 2017 г., максимальное количество сроков – 22032; холодный период ноябрь 2017 г. – март 2018 г., максимальное количество сроков – 21744

Table 8. Total number of forecast fields (times) with the objects identified after spatial smoothing with convolution radius of 9 grid pints using the isoline of 1 mm/h and with sizes of not less than 1600 grid points. Warm period is May-September 2017, the maximal number of times is 22032; cold period is November 2017 – March 2018, maximal number of times is 21744

 

ДМРЛ

 

Наблюдение

Прогноз (мин)

30

60

90

120

тепл

холод

тепл

холод

тепл

холод

тепл

холод

тепл

холод

RAKU

2843

1387

2550

1495

2561

1527

2553

1550

2561

1564

RATL

2378

847

2045

881

2045

908

2046

916

2018

931

RAVN

2147

130

*

*

*

*

*

*

*

*

RAVO

2777

506

2688

619

2728

655

2744

666

2748

687

RUDB

2246

721

1950

822

1936

858

1918

885

1932

884

RUDK

3432

379

3287

369

3266

402

3268

428

3237

445

RUDL

3240

691

3036

798

3036

744

3052

750

3088

753

RUDN

3520

626

3431

671

3456

694

3461

732

3427

767

RUWJ

3755

201

3533

235

3544

254

3554

266

3524

277

 

 

Из-за сбоев в информационных потоках в теплый период количество прогнозов оказалось примерно на 10 % меньше количества наблюдений [21], поэтому для корректного сопоставления следует вычесть из наблюдений эти 10 %: например, 2843 заменить на 2559 = 2843 – 284 и т.д.   

По столбцам "наблюдение" отчетливо видно, насколько разнообразно количество (fc) отбираемых полей для различных локаторов в оба периода: так, для теплого периода минимальное fc равно 2021, а максимальное – 3379; для холодного периода эти числа равны, соответственно, 201 и 1387 (исключая р/л Внуково). Естественно наследование такого различия в прогнозах. Выявленное разнообразие значений fc затрудняет обобщения по всем локаторам. Например, осреднение этих величин дает превышение значений fc в теплый период над аналогичными значениями в холодный период примерно в 4.5 раза, в то время как по локаторам это отношение колеблется от 2.1 до 9.1 в семи случаях и 18.7 – в одном (Валдай, RUWJ). При этом локатор Валдай обеспечивает максимальное значение fc по теплому периоду и минимальное значение в холодный период и в наблюдениях, и во всех прогностических сроках. Точно так же локатор Брянск (RUDB) обеспечивает минимальное значение fc в теплый период, а локатор Курск (RAKU) обеспечивает максимальное значение fc в холодный период.

При переходе от наблюдений к прогнозу 30 мин значение fc в теплый период меняется примерно на 100–200 полей: возрастает для пяти локаторов, а падает – для трех локаторов; в холодный период меняется примерно на 70–100 полей, при этом возрастая для семи и падая для одного локатора.

 

5. Учет пропусков в узлах полей и формирование выборок разных типов

В усовершенствованной системе наукастинга осадков, построенной на новом расчетном ядре pySTEPS, используются 1) расширенный набор форматов данных, 2) блок быстрого двумерного преобразования Фурье без ограничений на размер расчетного поля, 3) новая технология генерирования ансамбля и 4) новая методика обработки пропусков. Основные результаты испытаний системы в теплый период 2020 г. представлены в [22, 23]. Архив испытаний состоял из около 22000 ансамблевых прогнозов на 150 минут в 10-минутной дискретности, при этом каждый из 15 прогностических сроков был представлен ансамблем из 10 полей.

В отличие от предыдущих испытаний (2017–2019 гг.), в которых пропуски (пустые точки) замещались малыми положительными числами и проблема отсутствующих данных в принципе не возникала, необходимость обработки пропусков в 2020 г. привела к соответствующему изменению и дополнению стратегии верификации прогнозов с учетом сложившейся методологии обработки пропусков [19, 45, 49, 50]. Для оценки временного и пространственного распределения пустых точек проведен анализ полного архива наблюдений и полного архива ансамблевых прогнозов для одной заблаговременности (60 мин). Результаты анализа качества наукастинга, проведенного с помощью средств языка R, представлены в таблицах и графиках. Пространственный анализ производился с помощью квантильного анализа и соответствующих пространственных квантильных функций. 

 

5.1. Статистика пропусков в архивах наблюдений и прогнозов

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

Матрицу, полученную суммированием всех масок и которую назовем суммарной маской, проанализируем с помощью квантилей. Функция quantile языка R имеет на входе список порогов вероятностей, а на выходе набор квантилей, рассчитанных с помощью линейной интерполяции эмпирической функции распределения по методу Хайндмана [43]. В частности, последовательность из 11 порогов 0.0, 0.1, ..., 1.0 выводит такое же количество квантилей, или десятков процентов, от минимального значения (порогу 0.0 соответствует 0 %) до максимального (порогу 1.0 соответствует 100 %). Матрица, полученная произведением всех масок, состоит из единиц и нулей, при этом, если в точке поля имеется хотя бы одна дата с пропуском, результат будет равен нулю. Сумма единиц в такой матрице эквивалентна оценке площади (или размера) области без пропусков за весь период верификации.

Рассмотрим вначале суммарную маску: наибольшее значение в точке равно 22026, минимально возможное значение равно единице (этот факт, сам по себе довольно любопытный, тоже встречается). При наличии в полях в среднем более 70 % точек с числовыми значениями в пропорции к полному размеру, существуют поля с заметно меньшим количеством "непустых" точек (не более 50 %). Наибольшее количество точек поля с ненулевыми значениями равно 3193798, т. е. примерно в 76 % точек объединенного поля хотя бы раз встречается числовое значение.  

Выделим следующие особенности. Во-первых, более 95 % значений точек обеспечено не менее 18372 наблюдениями, а это примерно 95 % от 3193798 "непустых" точек. Одновременно с этим среди 5 % точек могут оказаться точки с чрезвычайно скудным обеспечением: менее 1000 числовых значений было обнаружено во временных рядах 2716 точек поля, менее 10 значений – в 1233 точках и ровно одно числовое значение – в 242 точках. Непросто определить, какое количество точек поля с "непустыми" значениями считать несущественным, и это зависит от формы области: связная область в 2500 точек может занимать квадрат со стороной 50 км, а протяженную линию – в 2500 км.

При шаге вероятности 0.05 квантили распределения значений в точках суммарной маски с ненулевыми значениями собраны в табл. 9.

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

Для удобства изображения больших чисел данные на обеих панелях рис. 7 приводятся в долях от максимального количества возможных случаев (22026). Видно четкое проявление кругов обзоров одинакового цвета для одного локатора (одинаковое заполнение числами всех точек поля, кроме теневых областей от местников), при этом наличие перекрывающихся кругов обзора приводит к подавлению теневых областей, кроме, возможно, не устраненных лучей и конусов р/л Внуково, Нижнего Новгорода, Краснодара (рис. 7а). На карте долей в последнем десятичном интервале (0.90, 0.92, ..., 0.98, 1.00) заметно отсутствие кругов р/л Воейково (0.8–0.9), Вологды (0.5–0.6), Белгорода (0.8–0.9).

 

Таблица 9. Квантили распределения количества числовых значений в 22026 объединенных полях радиолокационных наблюдений

Tabe 9. Distribution quantiles of the numerical values in 22026 merged radar fields

 

0%    5%      10%     15%     20%     25%     30%    35%     40%     45%    50%   

1    18347  21000  21468  21812  21855  21884  21911  21936  21986  21994

     ...продолжение...

  55%     60%     65%     70%     75%    80%    85%     90%     95%    100%

21998  22013  22018  22023  22023  22026  22026  22026  22026  22026

Примечание. Шаг частоты равен 0.05 (5 %); минимальное значение соответствует квантилю 0 %, что эквивалентно единственному наблюдению в некоторых точках поля; максимальное значение соответствует квантилю 100 %, что эквивалентно наличию в этих точках максимально возможного количества в 22026 наблюдений.

 

 

а)

б)

Рис. 7. Карты суммарных масок наличия числовых данных в круговых зонах обзора радиолокаторов усовершенствованной системы наукастинга осадков. Количество представлено в виде долей от максимально возможного значения 22026: весь диапазон значений (от 1 до 22026) (а); значения выше квантиля q(90%) (б).

Fig. 7. The maps of the summary masks of numerical data availability in the radar scanning circle for the improved nowcasting system. The number is represented as the fraction of the maximal possible number of 22026: the whole range from 1 to 22026 (а); values above quantile q(90%) (б).

 

Рис. 8. Поточечное произведение всех масок. В этих точках количество значений максимально (22026). Количество точек с единицами составляет 721083/3193798 = 23 % от количества всех "непустых" точек.

Fig. 8. Grid-by-grid product of all the masks. The number is maximal at these points (22026). The number of points with unity (1) is 721083/3193798 = 23 % of all "not-empty" points.

 

 

Можно сделать вывод, что обеспечение данными по всей области наукастинга вполне удовлетворительное, но возникают вопросы частного и общего характера. Например, можно ли считать обеспечение данными области обзора р/л Вологды неудовлетворительным при том, что доли от 0.5 до 0.6 означают суммарно не менее двух месяцев числовых наблюдений (144 срока в сутки)? И каким образом вообще можно оценить влияние асинхронных доступов в точках поля на пространственную согласованность показателей качества?

На карте произведения масок (рис. 8) отчетливо видно, что основная порция максимально обеспеченных точек приходится на области пересечения зон обзора локаторов в центральной части России (область ЦФО) и в районе локаторов Ставрополь и Сочи.

 

 

Пропуски в полях прогностических ансамблей. Имеющиеся в полях наблюдений пропуски естественно наследуются полями прогноза, при этом в основном размножаясь с ростом заблаговременности. Рассмотрим маски в ансамблях наукастов. В табл. 10 показаны сопоставленные ("спаренные") в прогнозе на 60 мин объединенные и нумерованные (от 1 до 21995) пары полей наблюдений (префиксы nc-файлов – "PRCP-ETR") и прогностических ансамблей (префиксы tiff-файлов – "steps_mvp") с указанием количества пропусков (первое значение в строке с номером и вторая строка ниже) и количества числовых значений (второе значение в нумерованной строке и третья строка ниже). 

 

Таблица 10. Количество пропусков и числовых значений в сопоставленных в прогнозе на 60 мин парах полей наблюдения (префиксы nc-файлов – PRCP-ETR) и прогностических ансамблей (префиксы tiff-файлов – steps_mvp)

Table 10. Number of missing data and numerical values in the observation fields (prefix PRCP_ETR) and forecast ensembles (prefixes of tiff-files step_mvp) pairs for the lead time of 60 min

Примечание. Количество пропусков – первое значение в строке с номером (по первой паре: 1005711) и вторая строка ниже (десять значений от 1018623 до 1031277 и среднее по ансамблю – 1023748.6). Количество числовых значений – второе значение в нумерованной строке (3190890) и третья строка ниже (от 3177978 до 3165324 и среднее по ансамблю – 3172852.4). Последний столбец – разности соответствующих значений количества пропусков в прогнозе и в наблюдении: 18037.6 = 1023748.6 – 1005711. Для контроля приводится аналогичная разность для количества числовых значений, которая в сумме с предыдущей разностью ожидаемо равна нулю.

 

 

Отметим самые общие особенности: 1) наблюдается большое разнообразие в количестве пропусков как в полях наблюдений, так и в полях ансамблей, включая поля внутри одного ансамбля, 2) как правило, количество пропусков в ансамблях превышает, иногда существенно, количество пропусков в соответственном поле наблюдения. Первая особенность полностью зависит от еще не проясненных условий формирования объединенных полей на основе полей отдельных локаторов. Вторая особенность представляется еще более интересной из-за значительного разброса разностей от нескольких сотен (например, 432 в табл. 3) до нескольких сотен тысяч (например, 281000 в той же таблице). Имеются случаи, когда количество пропусков в наблюдениях превышает их среднее или даже максимальное по ансамблю количество пропусков.

Так как при расчете показателей качества пустая точка любого из 11 полей (наблюдения и десяти полей ансамбля) исключается из результирующего поля показателя, то более содержательной оценкой различия следует считать разность между максимальным (а не средним) количеством пропусков в ансамбле и количеством пропусков в поле наблюдений. При этом, например, по первой паре вместо 18037 получим 1031277 – 1005711 = 25566.

Суммарная характеристика табл. 10 отражена в квартильных оценках в табл. 11. Максимальное количество пропусков в полях наблюдений (2098901) почти совпадает с минимальным количеством числовых значений (2097700), однако в среднем (по mean и median) количество пропусков занимает около четверти, а количество числовых значений – около трех четвертей объединенного поля (4196601). Максимальное количество пропусков в прогностических ансамблях изменяется от 144592, т. е. намного меньше минимального количества пропусков в соответственных полях наблюдений (1005542), до 2120232, т. е. больше максимального количества пропусков в полях наблюдений (209890). В более важных для расчета показателей согласованных разностей (строка mxnas-onas) размах оказывается существенно большим: от -1066062 до +1114475.

 

Таблица 11. Квартильные характеристики и среднее значение (mean) количества пропусков (nas), числовых значений (vld) в полях наблюдений (obs) и прогнозов на 60 мин (frc), максимальное значение в ансамблях (mx), а также разность между максимальным количеством пропусков в ансамбле и количество пропусков в "спаренном" поле наблюдений (mxnas-onas)

Table 11. Quartile characteristics and the mean of the missing data number (nas), numerical values (vld) in observation fields (obs) and 60 min lead time forecasts (frc), the maximal value in the ensembles (mx), and the difference between the maximum number of missing data in the ensemble and in the merged observation field (mxnas-onas)

 

    Min          1st Qu.      Median      Mean      3rd Qu.       Max

obs_nas

 1005542   1005783   1053784   1067155   1106053   2098901

obs_vld

 2097700   3090548   3142817   3129446   3190818   3191059

mx_frc_nas

   144592   1055335   1089212   1100554   1139304   2120232

mxnas-onas

-1066062       22095       32654       33399       45431   1114475

 

Как говорилось выше, анализ пропусков в согласованных полях прогнозов и наблюдений важен для оценки потерь значимых точек в результирующем поле показателя качества. Не менее интересен вопрос о "путешествии" в прогнозах тех пустых точек начальных полей, которые сохранились в модели STEPS после двумерной Фурье-фильтрации перед входом в прогностический блок. Для этого потребуется сравнивать количество пропусков в серии прогностических полей одной заблаговременности (например, 1 час) и количество пропусков в паре соответствующих начальных полей модели STEPS. Так как архив прогнозов по срокам не имеет разрывов во времени, то можно воспользоваться табл. 9 и сравнить прогностические поля с отстоящими на час раньше полями наблюдений. Статистика разностей между максимальным количеством пропусков в ансамблях и количеством пропусков в двух соответствующих начальных полях (a1, a2) собрана в табл. 12.   

 

Таблица 12. Разности между максимальным количеством пропусков в ансамбле прогнозов на 60 мин (mxnas) и количеством пропусков в соответствующих начальных полях (a1 и а2)

Table 12. Differences between the maximal number of missing data in the 60 min forecast ensemble (mxnas) and in the corresponding initial fields (a1 and a2)

 

    Min         1st Qu.     Median    Mean    3rd Qu.     Max

mxnas-а1nas

-1051771    23611      32508     33414    43320   1114585

mxnas-а2nas

-1051936    24036      32437     33410    42695   1074679

 

 

Данные табл.12 подобны данным последней строки табл.11. Отметим редкое исключение: количество одновременных отрицательных значений относительно первого и второго поля оказывается равным 103 – только в этих случаях (из порядка 20 тысяч прогнозов) максимальное количество пропусков в прогностических ансамблях оказалось меньше пропусков в обоих начальных полях.

 

5.2. Формирование выборок для верификации и подсчет количества допусков

Так как сведения о количестве числовых значений (допусков) важны для оценок статистической значимости рассчитываемых показателей, то необходимо учитывать не только и не столько "естественные" потери из-за особенностей наблюдательных систем и прогностической модели, но и изменение количества "допущенных" значений в зависимости от способов формирования этих выборок. Будем обозначать допуски через vld (valid); использованное в верификации среднее по ансамблю прогностическое поле будем обозначать через AVR; как правило, типы выборок будут помечаться прописными буквами в тексте (Е, ЕТ и т. д) и строчными в таблицах и рисунках (е, еt и т. д).

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

Верификация прогнозов проводилась на парах чисел прогноз-наблюдение (frc, obs), удовлетворявших условиям отбора допусков при предварительном замещении константы отсутствия NA отрицательным числом -999. Допустимость определяется на конъюнкциях (&) и дизъюнкциях (|) неравенств следующим образом:

1. (frc ≥ 0) & (obs ≥ 0) – тип "Е";

2. (frc > 0 & obs ≥ 0) | (frc ≥ 0 & obs > 0)  – тип "ET";

3. (frc > 0.05 & obs ≥ 0) | (frc ≥ 0 & obs > 0.05)  – тип "TT";

4. (frc > 0) & (obs > 0) – тип "T".

Условие (1) отбирает все числовые пары значений, включая нулевые ("без осадков"); при условии (2) хотя бы одно из чисел должно быть положительно ("осадки или в прогнозах, или в наблюдениях", "или" – не исключительное, т. е. "осадки и там, и там" входит в список "допустимой пары"); при условии (3) отбираются значения, не меньшие порога 0.05 мм/ч; условие (4) требует ненулевые значения в обоих числах ("осадки одновременно в прогнозах и в наблюдениях").

Во всех четырех случаях пропуски исключаются автоматически благодаря замещению числом -999. Количество допустимых пар подсчитывается на каждом шаге цикла и в конце цикла используется в операциях осреднения при расчетах средних оценок качества, а соответствующая матрица сумм сохраняется в виде поля исходного размера (nx, ny), nx×ny=4196601.   

Рассмотрим матрицы сумм числовых значений с этими четырьмя условиями в квартильном представлении, включая оценку среднего количества (mean) (табл. 13). Прокомментируем данную таблицу. Во-первых, видно, что среднее количество пропусков равно ~1 млн точек, количество числовых значений ~3 200 000 точек, и эти значения мало зависят от типов отбора. Во-вторых, изменение квартилей и средней в зависимости от заблаговременности можно считать незначительным, но заметно изменение квартилей при смене типов отбора допусков: например, для прогноза на 60 мин 75%-й квартиль меняется последовательно 21992 ("е"), 3943 ("еt"), 3014 ("tt") и 1942 ("t").

 

 

Таблица 13. Распределение в квартилях числовых значений по количеству "допустимых" точек поля с указанием количества пропусков (NAs) и числовых значений (vld) для различных заблаговременностей (lead в мин). Прогноз по среднему полю ансамбля (AVR)

Table 13. Quartile distribution of numerical values of valid points with the numbers of missing data (NAs) and valid numerical data (vld) for different lead times (lead in min). Ensemble average nowcast (AVR)

 

Квартили количества числовых значений
 в точке поля

 Пропуски   Допуски

lead

min   q25        med       mean        q75         max

          NAs         vld

тип отбора пар числовых значений "е"

030

060

090

120

150

 1     21766     21942     21412     21998     22000

 1     21733     21930     21188     21992     21995

 1     21556     21920     20893     21984     21990

 1     21288     21860     20529     21981     21987

 1     20848     21818     20095     21976     21984

1004617     3191984

1004764     3191837

1004919     3191682

1005011     3191590

1005096     3191505

тип отбора пар числовых значений "еt"

030

060

090

120

150

 1       1897       2631      2605      3226      5016

 1       2362       3275      3200      3943      5944

 1       2722       3760      3655      4504      6608

 1       3025       4155      4031      5022      7353

 1       3290       4500      4355      5534      8043

1004696     3191905 1004832     3191769

1004965     3191636

1005075     3191526

1005133     3191468

тип отбора пар числовых значений "tt"

030

060

090

120

150

 1       1510       21