Исследование Земли из Космоса, 2023, № 5, стр. 23-35

ВВЕДЕНИЕ

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

Спутниковые оценки количества осадков могут быть получены на основе измерений различных приборов. Методы в основном разделены на три главные категории, основанные на типе наблюдений, а именно: методы обработки данных, полученных в видимом, инфракрасном диапазонах; в микроволновом диапазоне излучения; а также методы, основанные на их комбинировании. В настоящей работе будут рассмотрены только микроволновые методы, а именно радиометрические наблюдения. Радиометрические сканеры/зондировщики, установленные на спутниках, приобрели наибольшую популярность при измерении интенсивности осадков в глобальном масштабе. Их основное преимущество перед локаторами заключается в малом энергопотреблении, что позволяет проводить измерения в течение многих лет. Как пример, наиболее известный прибор TMI (TRMM Microwave Imager) проработал 15 лет.

В настоящее время существует много радиометрических приборов, используемых для оценки интенсивности осадков. Это ATMS (Advanced Technology Microwave Sounder), MWRI-1 (Micro-Wave Radiation Imager-1), SSMIS (Special Sensor Microwave – Imager/Sounder), AMSR-2 (Advanced Microwave Scanning Radiomete) и система GMI-Core (Global Precipitation Measurement Microwave Imager) в которой собираются все доступные радиометрические наблюдения (Kummerow et al., 2015; Zabolotskikh, Chapron, 2015; Zhang et al., 2018; Surussavadee, Staelin, 2010). Среди отечественных, есть серия приборов МТВЗА (Модуль Температурного и Влажностного Зондирования Атмосферы) которые запускались на КА “Метеор-М” № 1 (2009 г.), “Метеор-М” № 2 (2014 г.) и “Метеор-М” № 2-2 (2019 г). В докладе (Бухаров и др., 2010) было упомянуто наличие у их авторов некоторых методик оценки интенсивности осадков над водной поверхностью. Однако, отсутствие в открытом доступе каких-либо публикаций по этому вопросу не позволяет судить о возможности их применения к длинным рядам данных МТВЗА-ГЯ (разных серий в том числе) и оценки численных показателей их качества, как, например, в (Ferraro and Marks, 1995). Таким образом, необходимо продолжать исследования и разработку методов восстановления интенсивности осадков по данным МТВЗА-ГЯ, что и является основной целью настоящей работы.

В микроволновой части спектра излучение поверхностью Земли изменяется наличием частиц, размером с дождевые капли. Для определения осадков используется два физических процесса: (1) излучение гидрометеоров, которое приводит к увеличению микроволнового излучения, и (2) рассеяние, что приводит к уменьшению интенсивности излучения. Величина этих эффектов зависит от размера и количества частиц, частоты принимаемого излучения (Kidd et al., 2010). Алгоритмы определения осадков, применительно к пассивным микроволновым наблюдениям, могут быть разделены на основе использования микроволновых каналов или методологически. Выделяют три класса алгоритмов: (1) класс “излучение” – этот тип алгоритмов использует “низкочастотные” каналы для определения увеличения излучения по сравнению с “холодным” излучение поверхности океана (Chang et al., 1999), (2) класс “рассеяние” – этот тип алгоритмов связывает осадки с уменьшением излучения, вызванным рассеянием на частицах (Ferraro and Marks, 1995), и (3) это тип “многоканальной инверсии” – комбинирование первых двух классов (Kummerow et al., 2001). По методологической классификации, алгоритмы могут быть отнесены к одной из двух групп: эмпирические методы – разрабатываются на основе экспериментальных данных и физические методы, которые стараются минимизировать разницу между модельными и экспериментальными данными.

В настоящей работе предложено начать с адаптации простой регрессионной модели восстановления интенсивности осадков над поверхностью океана, модель ALG'85 (Ferraro, 1997) используя только данные прибора МТВЗА-ГЯ. Это позволит исследовать возможность восстановления интенсивности осадков по измерениям МТВЗА-ГЯ и определить ряд проблем/задач, которые требуют дополнительных исследований.

МОДУЛЬ ТЕМПЕРАТУРНОГО И ВЛАЖНОСТНОГО ЗОНДИРОВАНИЯ АТМОСФЕРЫ (МТВЗА-ГЯ)

Приборы серии МТВЗА похожи на свои зарубежные аналоги. Основным отличием серии МТВЗА является необычный, для приборов с коническим сканированием, угол визировании (53.3°) и соответственно угол падения в 65°. Для большинства зарубежных аналогов угол падения составляет порядка 53°. Такое отличие приводит к необходимости разрабатывать новые (или переделывать уже имеющиеся) алгоритмы обработки данных, для того что бы конечный пользователь мог получить качественную информацию из нового независимого источника. (Ермаков и др., 2021).

Радиометрические каналы сканера МТВЗА-ГЯ включают рабочие частоты в окнах прозрачности атмосферы 10.6, 18.7, 23.8, 31.5, 36.7, 42, 48 и 91.65 ГГц, а зондировщика – в линиях поглощения кислорода 52–58 ГГц и водяного пара 183 ГГц. Подробное описание можно найти в работах (Болдырев и др., 2008; Чернявский и др., 2018).

Спутниковые данные МТВЗА-ГЯ поступают в ЦКП (Центр коллективного пользования системами архивации, обработки и анализа данных спутниковых наблюдений ИКИ РАН) “ИКИ-Мониторинг” в виде сформированных потоков данных и сопутствующей информации звездных датчиков (Ермаков и др., 2021). Результаты сохраняются в архивах ЦКП “ИКИ-Мониторинг” (Лупян и др., 2019; http://ckp.geosmis.ru/), что обеспечивает автоматическое формирование наборов данных для дальнейшей обработки и отображения в картографическом веб-интерфейсе системы Вега-Science (http://sci-vega.ru/).

ТЕХНИКА ВЫДЕЛЕНИЯ СВОБОДНОЙ ВОДНОЙ ПОВЕРХНОСТИ

Восстановление интенсивности осадков над поверхностью воды подразумевает наблюдение сцен свободной водной поверхности (без суши и льда). Для выделения суши была взята топографическая карта с шагом 60 угловых секунд. (NOAA National Centers for Environmental Information. 2022: ETOPO 2022 15 Arc-Second Global Relief Model. NOAA National Centers for Environmental Information. https://doi.org/10.25921/fd45-gt74. Accessed [15.01.2022]). По указанной топографической карте была создана маска для картирования данных МТВЗА-ГЯ по типам: суша, берег и вода. Так как пространственное разрешение самого низкочастотного канала МТВЗА-ГЯ (10.6 ГГц) составляет порядка 100 км по одной из осей, то отступ от берега был увеличен на 1 градус в сторону воды для того чтобы все анализируемые измерения находились над поверхностью воды.

Картирование ледовой области представляет определенные трудности, связанные с пространственным и временным изменением площади ледового покрытия в океане. Одним из важных источником информации для оперативных задач мониторинга ледяного покрова являются поля сплоченности льда (англ. sea ice concentration – SIC). Одно из недавних исследований по восстановлению SIC по данным МТВЗА-ГЯ представлено в работе Заболотских и др., 2022. Часть предложенного алгоритма использована для текущей задачи картирования ледовой области.

Традиционно, для восстановления SIC используют данные полученные на частотах вблизи 90 ГГц, благодаря чему достигается максимально возможное на сегодняшний день разрешение при использовании радиометрии. Однако, у прибора МТВЗА-ГЯ не работает канал с горизонтальной поляризацией на частоте 91.65 ГГц, поэтому использовать общепринятые алгоритмы картирования ледовой области не представляется возможным. В работе (Заболотских и др., 2022) предложен алгоритм, по которому возможно определить концентрацию ледяного покрова используя частоты 10.6 и 36.7 ГГц. В текущем исследовании нет необходимости определять концентрацию льда, а достаточно определить наличие льда на поверхности, поэтому будет использована только поляризационная разница (PD – разница между вертикально и горизонтально поляризованным излучением на одной частоте) на частоте в 10.6 ГГц. Использование указанной частоты ухудшает пространственное разрешение, однако, по оценкам (Заболотских и др., 2022) позволяет восстанавливать SIC с меньшими погрешностями, а, следовательно, и лучше будет проводиться картирование ледовой области. Указано, что при PD < 120 К (10.6 ГГц) начинает появляться лед, поэтому значение в 120 К выбрано как пороговое.

Известно, что максимальная площадь ледового покрова Арктики обычно наблюдается в марте, а Антарктики в конце сентября. Для этих периодов были проанализированы данные полученные из системы Remote Sensing Systems DMSP SSM/I or SSMIS (Wentz et al., 2012). По полученные картам была построена маска максимальной площади ледового покрытия регионов Арктики и Антарктики. К построенной маске максимальной площади ледового покрытия добавлен запас в 2 градуса по широте. Если наблюдения попадают в область, где возможен лед, то проходят проверку поляризационной разницы, если PD < 120 К, то данное измерение маркируется как лед и не используется при анализе осадков.

ИНДЕКС РАССЕЯНИЯ. ПОДБОР УРАВНЕНИЯ

Для определения интенсивности осадков используется два четких процесса: (1) излучение дождевых капель, которое приводит к увеличению микроволнового излучения на частотах ниже 40 ГГц, и (2) рассеяние, вызванное осадками, что приводит к уменьшению интенсивности проходящего через осадки излучения на частотах выше 80 ГГц. Значения указанных пороговых частот радиотеплового излучения условное, и связно с работой радиометрических систем, используемых для мониторинга поверхности, в окнах прозрачности атмосферы. Величина эффектов излучения и рассеяния на дождевых каплях зависит от размера и количества частиц. Однако возможности измерения осадков определяются также и свойствами подстилающей поверхности. В работе (Grody, 1991) проведено обширное исследование классификации снежного покрова и осадков на основе измерений прибора SSMI (Special Sensor Microwave Imager), и показано, что над водной поверхностью классифицировать осадки проще, чем над поверхностью суши. Это связано с неоднородностью излучающей способности суши по сравнению с практически однородной излучающей способностью морской воды.

В работе (Grody, 1991) используется Индекс Рассеяния (англ. Scattering Index) SI как положительная разница между измерениями на частотах 22 и 85 ГГц. Данный индекс характеризует наличие веществ в атмосфере и на поверхности, которые рассеивают восходящее радиотепловое излучение. Так же было отмечено, что для более точной идентификации рассеяния необходимо привлекать измерения на других частотных каналах и что оценки должны быть выполнены на большом наборе данных, включающем в себя различные состояния поверхности и атмосферы.

Применительно к радиометрическим каналам прибора МТВЗА-ГЯ индекс рассеяния для частоты 91.65 ГГц вертикальной поляризации (V), далее просто SI, может быть записан как

где F – это функция, выраженная в виде суммы с весовыми коэффициентами радиояркостных температур на частотах отличных от 91.65 ГГц, и характеризующая каким могло бы быть излучение на частоте 91.65 ГГц в отсутствии рассеяния. $T_{{91.65}}^{V}$ – радиояркостная температура на частоте 91.65 ГГц вертикальная поляризация.

У прибора МТВЗА-ГЯ всего доступно 24 радиометрических канала. Однако, на некоторых канала есть сбои и высокий уровень шумов, они не используются (это каналы $T_{{91.6}}^{H},$ $T_{{42}}^{{V,H}},$ $T_{{48}}^{{V,H}}$). Так же не используются каналы в линии поглощения кислорода. Поэтому для анализа индекса рассеяния SI будет использовано всего 11 каналов данных. Это каналы $T_{{10.6}}^{{V,{\kern 1pt} H}},$ $T_{{18.7}}^{{V,{\kern 1pt} H}},$ $T_{{23.8}}^{{V,{\kern 1pt} H}},$ $T_{{31.5}}^{{V,{\kern 1pt} H}},$ $T_{{36.7}}^{{V,{\kern 1pt} H}}$ и $T_{{91.65}}^{V}.$ Для разработки соотношения для функции F используются все данные за 2020 г., усредненные за каждый месяц. Усреднение было сделано для уменьшения объема точек и многократного сокращения вычислительных ресурсов. Усреднение проводилось для всех измерений над открытой поверхностью океана свободной ото льда в диапазоне ±0.5° по широте для всех значений долгот. Таким образом, за каждый месяц получены средние значения радиояркосных температур для каждого радиометрического канала МТВЗА-ГЯ в зависимости от широты.

Для составления регрессионного соотношения функции F изначально были использованы все радиометрические каналы данных в “окнах прозрачности” атмосферы.

где: а_i – коэффициенты регрессии. Под радиометрическим каналом Т_i понимается измерение на одной частоте и на одной поляризации. Коэффициенты регрессии вычисляются методом наименьших квадратов по методике, указанной в статье (Сазонов и др., 2020). Общий коэффициент корреляции между радиояркостной температурой $T_{{91.65}}^{V}$ и функцией F составил 0.9961, а СКО не более 2 К.

Соотношение (2) адекватно связывает радиояркостную температуру $T_{{91.65}}^{V}$ с ее расчетом на основе других радиометрических каналов прибора МТВЗА-ГЯ, однако, избыточно и может быть упрощено с помощью проверки коэффициентов регрессии по уровню значимости. Это позволит убрать из уравнения (2) факторы, влияние которых незначительно, либо полностью отсутствует, что упростит вид регрессионного соотношения.

Методика проверки коэффициентов линейной регрессии на статистическую значимость приведена в приложении статьи (Сазонов и др., 2020). Для коэффициентов регрессии (2) были рассчитаны значения t-критерия (t_a). Из таблицы критических точек распределения Стьюдента выбирается пороговое значение t(k, a = 0.01) = 2.58 (а – уровень значимости, k – число степеней свободы, при большом количестве экспериментальных данных, k принимается равным бесконечности). Таким образом, если t_a < 2.58 для какого либо коэффициента а_i, то этот коэффициент статистически незначим и может быть убраны из уравнения регрессии без потери качества модели. Результирующее уравнение будет иметь следующий вид:

где коэффициенты ${{а}_{0}}{\kern 1pt} - {\kern 1pt} {{a}_{8}}$ приведены в табл. 1.

На рис. 1 приведена диаграмма рассеяния усредненных за месяц радиояркостных температур на канале $T_{{91.65}}^{V}$ и их модельных по регрессии (3) значений. Корреляция составляет 0.9932. На рис. 2 приведены расчетные значения СКО, накопленные в диапазоне ±0.5 К. СКО, в среднем, составляет порядка 1.5 К. Большой значение СКО в области с температурами менее 235 К связан с малым количеством анализируемых данных в этой области.

На рис. 3 приведены измеренные значения радиояркостной температуры на канале $T_{{91.65}}^{V},$ ее модельный расчет на основе уравнения (3) и их разница на основе уравнения (1) – индекс рассеяния SI на частоте 91.65 ГГц. Индекс рассеяния используется для анализа и восстановления интенсивности осадков.

СОВМЕЩЕНИЕ ДАННЫХ РЕАНАЛИЗА И МТВЗА. ВЫБОР ОБЛАСТЕЙ С ОСАДКАМИ

Для исследования интенсивности осадков необходим большой набор наблюдений разных по интенсивности зон осадков в разных частях Земли. Одним из таких наборов является GPM IMERG (Integrated Multi-satellitE Retrievals for Global Precipitation Measurements) (Huffman et al., 2019). Указанный набор состоит из оценок осадков с различных пассивных микроволновых радиометров входящих в созвездие GPM (https://gpm.nasa.gov/missions/GPM/constellation). Осадки рассчитываются с использованием алгоритма Goddard Profiling Algorithm версии 2017 г. (GPROF2017) (Randel et al., 2020), затем приводятся к координатной сетке и проходят интеркалибровку по алгоритму Combined Ku Radar-Radiometer Algorithm (CORRA). Итоговый продукт – получасовые поля размером 0.1° × 0.1° (примерно 10 × 10 км).

Данные МТВЗА-ГЯ состоят из файлов записей восходящих и нисходящих полувитков. По отдельности работать с этими файлами достаточно сложно. Поэтому полученные данные были собраны в один массив за 1 день с разделением на восходящие и нисходящие витки. Совмещение набора GPM IMERG с данными МТВЗА-ГЯ проведено по времени, то есть, для каждого измерения МТВЗА-ГЯ выбирались данные из набора GPM IMERG с разницей по времени не более чем на ±15 мин. Таким образом, получены наборы данных об осадках в координатах широты, долготы и времени, соответствующих прибору МТВЗА-ГЯ. Затем, для удобства обработки и визуализации, полученные массивы всех данных накладываются на регулярную сетку 0.25° × 0.25° (рис. 3, карта интенсивности осадков). Собраны все имеющиеся данные МТВЗА-ГЯ за 2020 г. с 1 февраля по 31 декабря.

Выделение областей с осадками/индексом рассеяния производится в следующей последовательности: (1) На рассчитанные карты с индексом рассеяния SI накладывается маска при условии, что SI > 10. (2) Полученная маска преобразуется в контур, и выбираются замкнутые контуры с количеством точек периметра не менее 80 (данное значение получено опытным путем, и позволяет выделить средние и большие области с увеличенным рассеянием). (3) Определяются границы областей и к ним добавляется запас величиной в 1°. Полученные прямоугольные области используются для анализа связи индекса рассеянии и интенсивности осадков. За 2020 г. накоплено 3448 областей для восходящих витков и 4166 областей для нисходящих витков.

ПОСТРОЕНИЕ ЗАВИСИМОСТИ ИНТЕНСИВНОСТИ ОСАДКОВ ОТ ИНДЕКСА РАССЕЯНИЯ

На рис. 4 представлена анализируемая область с осадками, полученная по данным GPM IMERG, область с рассеянием, полученным по уравнению (1). Временное расхождение при совмещении данных составляет –1 мин, т.е. измерения, выполненные прибором МТВЗА-ГЯ проведены на 1 мин раньше, чем оценки интенсивности осадков по данным GPM IMERG. Как видно из представленных карт есть разница в расположении области рассеяния и области осадков. Это может быть связно с неточностью геопривязки наблюдений МТВЗА-ГЯ, о которой изложено в работе (Sadovsky, Sazonov, 2022). Также возможны ошибки сведения лучей диаграммы направленности МТВЗА-ГЯ для разных радиометрических каналов (этот вопрос еще не исследован). Таким образом, нет возможности наложить напрямую карты с рассеянием и осадками друг на друга.

Для получения зависимости интенсивности осадков от индекса рассеяния SI по выбранной области все значения, как осадков, так и рассеяния были построены в ряд и отсортированы по убыванию. Полученная зависимость показана рис. 4, и она соответствует результатам аналогичных исследований, показанных, например, в работе (Ferraro, 1997).

На рис. 5 приведены зависимости интенсивности осадков от индекса рассеяния SI. Статистика накоплена по всем выделенным областям для восходящих и нисходящих витков. Временное совмещение составляет ±1 мин. Как видно из рис. 5, результаты для восходящих и нисходящих витков в диапазоне осадков до 25 мм/ч практически идентичны. В диапазоне осадков более 25 мм/ч (и для SI > 50) данных очень мало и поэтому статистика не достоверна. Таким образом, для построения функциональной связи SI с интенсивностью осадков данные с обоих типов витков объединены.

Функциональную связь между индексом рассеяния и интенсивностью осадков I наиболее часто описывают степенным уравнением или полиномом. После исследования степенных, показательных и экспоненциальной функций в качестве аппроксимационных, было принято решение использовать полином 4-ой степени (4), который обладает наилучшим коэффициентом корреляции. Аппроксимация построена в предположении, что погрешности реанализа и ошибки моделирования SI отсутствуют. На рис. 6 приведены средние значения интенсивности осадков от SI, диапазон накопления Δ_SI = 2, и СКО, которое обуславливает 68% всех точек.

где значения коэффициентов полинома приведены в табл. 3 .

На рис. 7 приведено среднеквадратичное отклонение интенсивности осадков от модельных значений функции (4). СКО составляет в среднем 40–50% от восстанавливаемой величины.

ОБСУЖДЕНИЕ РЕЗУЛЬТАТОВ

На рис. 8 приведены три области с восстановленными интенсивностями осадков и соответствующие им значения осадков по данным GPM IMERG. Качественное сравнение указывает на существенную разницу между восстановленными осадками и данными GPM IMERG и на наличие смещения области осадков вниз для восходящих витков и вверх для нисходящих (здесь нисходящие витки не показаны). Указанная разница может быть связна с временным расхождением между наборами данных, которое для показанных областей составляет –1, 7 и 15 мин (соответственно, сверху вниз). Тогда на верхних областях мы должны наблюдать практически одинаковую картину распределения осадков. В принципе, так и есть если отбросить смещение и учесть погрешности из рис. 6 и 7.

Однако, большое значение СКО не может объяснить отмеченного смещения области восстановленной интенсивности осадков. Это можно объяснить, используя карты для индекса рассеяния SI на рис. 9. Над областью с большим рассеянием для восходящих витков есть темные участки, в которых значения нулевые и отрицательны. Для нисходящих витков эти участки находятся под областью осадков. В обоих случаях участки с отрицательными значениями SI повторяют контур области осадков. Исходя из формул (1) и (3), индекс рассеяния не может быть меньше нуля, однако, учитывая СКО для регрессии 3 (рис. 2) отрицательные значения возможны, но не более 5 К (в среднем), а получаемые значения SI могут достигать и –20 К, что физически не возможно.

Объяснить отрицательные значения SI можно только тем, что измерения для частоты 91.65 ГГц выполнены в другой точке пространства, нежели измерения для каналов с “низкой” частотой (10.6, 23.8 и 31.5 ГГц). Например, измерения на каналах с “низкой” частотой выполнены в области осадков. По формуле (3) проводится расчет F – модельной оценки излучения на частоте 91.65 ГГц. (Стоит учесть, что “низкочастотные” каналы тоже чувствительные к осадкам, хотя и в значительно меньшей степени, чем на частоте 91.65 ГГц, что можно отметить по 2-ому сверху графику на модельном расчете F (рис. 3)). Если луч диаграммы направленности канала 91.65 ГГц смотрит в туже точку, что и другие, тогда в зоне с осадками на данной частоте должно быть большое рассеяние (низкая радиояркостная температура) и результат по формуле (3) получится положительный. Если же диаграммы направленности канала 91.65 ГГц смотрит в другую, точку пространства, где осадков нет, тогда нет рассеяния на частицах осадков (высокая радиояркостная температура) и результат по формуле (3) получится отрицательный. В обратном случае, когда на канале 91.65 ГГц регистрируются осадки, а на других нет, будет положительное значение индекса рассеяния SI.

Предложенное объяснение, по-видимому, и определяет наблюдаемый сдвиг области рассеяния на рис. 9 и, как следствие, сдвиг области осадков на всех левых областях рис. 8. Таким образом, можно сделать вывод, что луч диаграммы направленности канала 91.65 ГГц захватывает другой участок пространства, нежели лучи для каналов (10.6, 23.8 и 31.5 ГГц). Более того, наблюдаемое смещение областей с осадками носит систематический характер. Все, вышесказанное, позволяет предположить, что лучи диаграмм направленности для различных частотных каналов прибора МТВЗА-ГЯ сведены некорректно.

ЗАКЛЮЧЕНИЕ

Полученные в работе результаты показывают, что разброс СКО достигает 50%, а коэффициент корреляции не превышает 0.75, по сравнению с результатами работы (Zabolotskikh, Chapron, 2015), где по данным прибора AMSR-2 интенсивность осадков восстанавливается с коэффициентом корреляции 0.8 и СКО = 1 мм/ч. Такие показатели не обеспечивают должную надежность связи радиоизлучения системы “океан–атмосфера” и интенсивности осадков, восстановленных по данным МТВЗА-ГЯ.

Конечно, интенсивность осадков, возможно, восстанавливать по радиометрическим измерениям МТВЗА-ГЯ. Однако на данном этапе работ точность неудовлетворительная, по сравнению с другими приборами и алгоритмами. Необходимо провести углубленный анализ параметров работы прибора – его геопривязку по всем каналам, сведение лучей диаграмм направленности, приведение пространственного разрешения к единому масштабу. И только после решения указанных вопросов вернуться к восстановлению интенсивности осадков и других геофизических параметров.