Методика построения климатической модели индекса атмосферной рефракции по данным радиозатменного зондирования Текст научной статьи по специальности «Математика»

Всероссийская открытая научная конференция «Современные проблемы дистанционного зондирования, радиолокации, распространения и дифракции волн» - Муром 2022

УДК 551.581 DOI: 10.24412/2304-0297-2022-1-282-289

Методика построения климатической модели индекса атмосферной рефракции по данным радиозатменного зондирования

А. В. Шмаков1, М. Е. Горбунов 1,2

1 Институт физики атмосферы им. А. М. Обухова РАН,

119017, Москва, Пыжевский пер., 3

E-mail: ldr@ifaran.ru

2 Гидрометцентр России,

123242, Москва, Б. Предтеченский пер., 11-13

E mail: gorbunov@ifaran. ru

В настоящее время радиопросвечивание атмосферы Земли с помощью навигационных спутников GPS широко применяется для восстановления параметров нейтральной атмосферы на высотах до 30-35 км. Накопленные данные за период более 10 лет проведения радиозатменного эксперимента COSMIC предоставляют важную информацию для понимания атмосферной динамики и изменения климата Земли. В этой работе мы разрабатываем методику построения климатической модели нейтрального атмосферного индекса рефракции. Представлена регрессионная модель зависимости индекса рефракции от высоты, широты, долготы и дня года. На основе пробных расчетов оптимизировали параметры модели. Создан пакет программ для массовой обработки восстановленных профилей индекса рефракции, на основе данных эксперимента COSMIC за период с 2007 года по 2016 год. Ключевые слова: радиозатменное зондирование, индекс рефракции, климатическая модель

Method for the development of a climatological model of the atmospheric refractive index from radio occultation data

A.V. Shmakov1, M.E. Gorbunov 1'2

1 A.M. Obukhov Institute of Atmospheric Physics RAS,

2 Hydrometeorological Center of Russia,

GPS radio occultation measurements are now widely used to retrieve the neutral atmospheric parameters in the altitude range up to 30-35 km. Accumulated data over a period of more than 10 years of the COSMIC radio occultation experiment provides important information for understanding atmospheric dynamics and climate change on Earth. In this work, we propose a method for the development of a climatological model of the neutral atmospheric refraction index. A regression model of the dependence of the refraction index on altitude, latitude, longitude and day of the year is presented. Based on trial calculations, the model parameters were optimized. We implemented a software package for bulk processing of retrieved profiles of refractive index, from the COSMIC radio occultation data for the period from 2007 to 2016. Keywords: radio occultation, refraction index, climatological model

ВВЕДЕНИЕ.

Радиозатменное зондирование заключается в восстановлении атмосферных параметров на основе измерений радиосигналов, прошедших через атмосферу Земли [1, 2]. В качестве источников сигнала используются высокоорбитальные (высота около 20000 км) навигационные спутники, такие как GPS (Global Positioning System) или ГЛОНАС (Глобальная навигационная спутниковая система). Радиосигналы проходят через ионосферу и нейтральную атмосферу и принимаются приемником, расположенным на низкоорбитальном спутнике (высота 700-800 км).

На первом этапе восстановления метеопараметров атмосферы из измеренных амплитуды и фазовой задержки сигнала определяют профиль угла атмосферной рефракции [3], который включает в себя ионосферную и нейтральную компоненты. Ионосферная коррекции, т.е. удаление ионосферной составляющей, позволяет восстановить профиль угла рефракции нейтральной атмосферы. Для этой цели используется зависимость ионосферного показателя преломления от частоты. Система спутников GPS излучает радиоволны на двух каналах (L1= 1.57542 ГГц и L2= 1.22760 ГГц). Для данного диапазона частот зависимость показателя преломления нейтральной атмосферы от частоты пренебрежимо мала, а показатель преломления ионосферы обратно пропорционален квадрату частоты. Используя линейную комбинацию углов рефракции из каждого канала при одинаковом прицельном параметре [4,] и его частоту, вычисляется профиль угла рефракции с удаленной ионосферной компонентой. Неточность линейного приближения, предложенного в [4], и влияния мелкомасштабной (менее 1 км) ионосферной турбулентности приводят к тому, что полученный профиль угла рефракции содержит остаточную погрешность ионосферной коррекции.

Эффективный метод подавления остаточных ионосферных шумов основан на применении статистической регуляризации [5]. Суть этого метода заключается в оценке наиболее вероятного профиля угла рефракции, полученного комбинацией среднестатистического модельного (априорного) и наблюдаемого профилей угла рефракции

В качестве априорной информации использовались модели CIRA (Committee on Space Research International Reference Atmosphere) [6] или MSIS (Mass Spectrometer and Incoherent Scatter Radar) [7]. Сравнение восстановленных параметров нейтральной атмосферы из радиозатменных данных с использованием этих моделей и данных Европейского Центра среднесрочного прогноза погоды (ECMWF) показали значительные расхождения на больших высотах [8].

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

В Датском Метеорологическом Институте этот подход был реализован при создании модели BAROCLIM (Bending Angle Radio Occultation Climatology) [11] на основе данных эксперимента COSMIC [12] за 2006-2012 годы. На конференции OPAC-2013 во время обсуждения предварительных результатов, полученных с помощью этой модели, были выявлены два недостатка [13]: на высотах выше 55 км для статистической обработки данных лучше использовать не средние профили угла рефракции, а медианы или робастные методы статистики и использование модели MSIS.

Учитывая недостатки модели BAROCLIM, в 2016 году авторами была разработана среднестатистическая модель угла рефракции BA-IAP (Bending Angle - Institute of Atmospheric Physics), построенная на основе обработки массива радиозатменных данных COSMIC за 2006-2013 годы [14]. Сравнение профилей угла рефракции показало, что модель BA-IAP лучше описывает среднестатистические углы рефракции по сравнению с моделью MSIS. Использование BA-IAP для восстановления профилей

показателя преломления показало, что выше 30-35 км для происходит систематическое уменьшение отклонения от данных ECMWF в сравнении с данными, восстановленными с использованием MSIS.

Целью данной работы является разработка методики построения климатической модели индекса рефракции на основе данных эксперимента COSMIC за 2007 - 2016 годы. Для этого будет использованы профили индекса рефракции, восстановленные из радиозатменных данных с помощью программы OCC, разработанной одним из авторов данной работы [15, 16]. Будет разработана регрессионная модель, где в качестве зависимой переменной будет профиль индекса рефракции, а предикторами будут зависимость от высоты, широты, долготы и дней года. Мы определим оптимальный состав предикторов и проведем пробный численный расчет коэффициентов модели. Будет подготовлен пакет программ для обработки массива данных эксперимента COSMIC за 10 лет.

Выбор предикторов и построение адаптивных функций.

Мы выбрали четыре основных предиктора для модели: высота к, широта р, долгота Я и день года /. Зависимость от высоты мы приближаем многочленами Чебышева Гп (гк), где п это порядок многочлена, и для этого делаем замену переменных

гк(к) = 2• (к~ко) -1, (1)

(км -к)

где км - максимальная высота наблюдения;

к - минимальная высота; п =10.

Широту и долготу мы берем в форме тригонометрических функций

1, оо8(р), $т(р), ооб(2 • р), Бт(2 • р), ооб(3 • р), Бт(3 • р); 1, ооб(Я), Бт(Я), ООБ(2 • Я), Бт(2 • Я) .

Зависимость от дня года берется в виде

1, г;

где г = 2 • -1. (2)

364

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

Полученное отношение для каждого предиктора анализировались и минимальные отбрасывались. Для примера на рис.1 представлен график значимости первых 20 полиномов Чебышева в качестве предикторов для зависимости индекса рефракции от высоты. Расчет делался для более чем 1000 восстановленных профилей индекса рефракции.

Адаптивные функции строятся как перекрестные произведения этих четырех наборов предикторов и в результате общее количество получается у = 10*7*5*2 = 700. Внутри каждой группы предикторов перекрестного перемножения не делается, так как полученные адаптивные функции не будут линейно независимыми.

...1 .... -__ МП 1111 НИ ______ .... INI ни im l_i_1 IUI Uli L...J _ IUI IUI ____ Uli 1_J _

О 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25

Рис. 1. График значимости предикторов в зависимости индекса рефракции от высоты. По оси абсцисс номер предиктора.

Описание модели.

Для построения климатической модели индекса рефракции мы используем профили, восстановленные из радиозатменных данных с помощью программы OCC на основе данных эксперимента COSMIC за 2007 - 2016 годы.

В общем случае задача регрессионного анализа заключается в нахождении взаимосвязи между зависимой переменной y и несколькими независимыми

переменными x (предикторами). Для этой цели строиться регрессионная модель в

которой мы аппроксимируем наблюдаемые величины y как линейную комбинация

адаптивных функций а1 предикторов x, где нижний индекс i определяет количество

наблюдений, а верхний индекс j количество адаптивных функций

yii«1 (xj (3)

j

обозначая

Kj =а (хг), (4)

запишем (3) в матричном представлении

У = K-a, (5)

где а1 коэффициенты регрессии, минимизирующие соотношение

(y - K - a)T (y - K - a) = min . (6)

В нашем случае количество наблюдений i во много раз больше числа адаптивных функций j поэтому система (3) сильно переопределена и решением (5) методом наименьших квадратов в матричной форме будет [ 17]

a = (KT • K)-1 (KT • y) . (7)

В результате перемножения KT • K получиться квадратная и симметричная матрица, обозначим ее A

Л, =S KP ■ KPj. (8)

p

Вектор, полученный в результате перемножения KT • y, обозначим z

= S Кл • y,. (9)

i

Запишем уравнение (7) для вычисления коэффициентов регрессии в виде:

a = A-1 • z . (10)

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

За десять лет проведения эксперимента COSMIC накопилось порядка N = 1010 профилей наблюдений, поэтому на практике для удобства и устойчивости машинного счета мы усредним матрицу B и вектор z на величину наблюдений N,

Л = ^S KP, • KPj = ^S>'(Хр)V(Xp), (11)

=1 s Kir y=1 (x)' *. (12)

Кроме этого, у нас 700 адаптивных функций, что потребует особого внимания и аккуратности при численных расчетах с такой матрицей и выбора подходящего алгоритма для ее обращения.

Один из самых удобных методов вычисления коэффициентов для решения задач методом наименьших квадратов основан на матричном разложении [18, 19], известном как разложение по сингулярным числам (Singular Value Decomposition, SVD) следующего вида

B = UAVT, (13)

где B любая вещественная матрица, U и V вещественные ортогональные матрицы такие, что диагональные элементы матрицы А имеют вид

У ^2 > ••• ^Уг >Уг+1 = ••• =Гп = 0, (14)

а не диагональные элементы равны 0, г - ранг матрицы В, п - ее размер. Если В не вырождена, то

У ~У2 ^••• ^Уп = 0

БУО-разложение является устойчивым, т.е. малым возмущениям матрицы В соответствуют малые возмущения матрицы Л и наоборот. Разработанные эффективные и устойчивые алгоритмы расчета БУО-разложения [20] позволяют использовать его для расчета обратной матрицы А-1. Подставляя А-1 в (10), мы получаем вектор коэффициентов а, что позволяет применять регрессионную модель (3) для получения профиля индекса рефракции в любой точке земного шара.

В нашем случае зависимость профиля индекса рефракции представляем в виде:

N (х) = ехр(^ а3 а (х)) + ехр(^ а3 а3 (х)) • ^ а3 (х) • 5а3 (15)

Решая его относительно да3 получим

£а3 (х) -да3 = ^ ^ -1, (16)

3 ехр(^ а3 а (х))

где N(x) -профиль индекса рефракции в атмосфере; х - набор предикторов; а3(х) - адаптивные функции;

а3 - искомые коэффициенты при адаптивных функциях; 3 - количество адаптивных функций.

Решение этого уравнения для да3 ищется итерациями. Для начального приближения предполагается, что ехр(^ а3пю3 (х)) просто экспонента. Вычислив да/ ,

I

на следующем шаге:

а(=а]0+да( . (17)

Подставляем полученное а3 в (16), вычисляем да3 и т.д., пока не выполниться соотношение (6).

Нами был проведен пробный численный расчет по данным за 540 дней. На рис. 2 и 3 представлены два произвольных восстановленных высотных профиля индекса рефракции (жирный черный цвет), полученные из модели профили (жирный красный цвет) и итерационные профили, получаемые в процессе счета.

Рис. 2. Высотные зависимости восстановленного (черный) и модельного (красный)

профиля индекса рефракции

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

Рис. 3. Высотные зависимости восстановленного (черный) и модельного (красный)

профиля индекса рефракции

На основе восстановленных профилей рефракции из данных эксперимента COSMIC была разработана регрессионная модель с предикторами высота, широта, долгота и дни года. Был определен оптимальный состав предикторов, для которых проведен пробный численный расчет коэффициентов модели. Сравнение модельных и восстановленных профилей индекса рефракции показало эффективность и адекватность построенной модели. Разработан алгоритм и подготовлен пакет программ для обработки массива данных эксперимента COSMIC за 10 лет.

Работа выполнена при поддержке гранта РФФИ №20-05-00189.

Литература

1. Kursinski E. R., Hajj G. A., Bertiger W. I., Leroy S. S., Meehan T. K., Romans L. J., Schofield J. T., McCleese D. J., Melbourne W. G., Thornton C. L., Yunck T. P., Eyre J. R., Nagatani R. N. Initial Results of Radio Occultation Observations of Earth's Atmosphere Using the Global Positioning System // Science, 1996, V. 271(5252), P. 1107-1110, doi:10.1126/science.271.5252.1107.

2. Kursinski E.R., Leroy S.S., Herman B. The GPS Radio Occultation Technique // Terr. Atmos. Oceanic Sci., 2000, V. 11. № 1. P. 53-114.

3. Rocken, C., Anthes, R., Exner, M., Hunt, D., Sokolovskiy, S., Ware, R., Gorbunov, M., Schreiner, W., Feng, D., Herman, B., Kuo, Y.-H. and Zuo, X. Analysis and validation of GPS/MET data in the neutral atmosphere // J. Geophys. Res., 1997,1V. 02, P. 29849-29866, doi :10.1029/97JD02400,

80

Индекс рефракции, N-ед.

Выводы.

4. Воробьев В. В., Красильникова Т. Г. Оценка точности восстановления атмосферного показателя преломления по измерениям доплеровского сдвига частоты на частотах, используемых в системе NAVSTAR // Известия АН СССР, Физика атмосферы и океана, 1993, Т. 29. № 5. C. 626-632.

5. Турчин В.Ф., Нозик В.З. Статистическая регуляризация решения некорректных задач // Изв.АН СССР. Физика атмосферы и океана. 1969. Т. 5, № 1. С. 14-18.

6. Fleming, E. L., Chandra, S., Barnett, J. J., and Corney, M. Zonal mean temperature, pressure, zonal wind, and geopotential height as functions of latitude // COSPAR International Reference Atmosphere 1986, 1990, Part II: Middle Atmosphere Models, Adv. Space Res., V. 10, P. 11-59.

7. Hedin, A. E.: Extension of the MSIS thermosphere model into the middle and lower atmosphere // J. Geophys. Res., 1991, V. 96, P. 1159-1172, doi:10.1029/90JA02125.

8. Gorbunov M. E., Shmakov A. V., Leroy S. S., Lauritsen K. B., COSMIC radio occultation processing: Cross-center comparison and validation // Journal of Atmospheric and Oceanic Technology, 2011, V. 28, No. 6, 737-751, doi: 10.1175/2011JTECHA1489.1

9. Ao, C. O., Mannucci, A. J., and Kursinski, E. R. Improving GPS radio occultation stratospheric refractivity retrievals for climate benchmarking // Geophys. Res. Lett., 2012, V. 39, L12701, doi:10.1029/2012GL051720.

10. Gleisner, H., Healy S. A simplified approach for generating GNSS radio occultation refractivity climatologies // Atmos. Meas. Tech., 2013, V. 6, P. 121-129, doi:10.5194/amt-6-121-2013.

11. Scherllin-Pirscher B., Syndergaard S., Foelsche U., Lauritsen K. B., Generation of a Bending Angle Radio Occultation Climatology (BAROCLIM) and its use in radio occultation retrievals // Atmos. Meas. Tech. Discuss., 2014, V. 7, P. 8193-8231, doi:10.5194/amtd-7-8193-2014.

12. Rocken, C., Y.-H. Kuo, W. S. Schreiner, D. Hunt, S. Sokolovskiy, and C. McCormick COSMIC System Description // Terr. Atmos. Oceanic Sci., 2000, V. 11(1), P. 157-186.

13. Gleisner H., Healy S., Refractivity and Temperature Climatology Retrieval Using Simplified Upper-Level Initialization, International Workshop on Occultations for Probing Atmosphere and Climate // 5-11 September 2013, Seggau Castle, Leibnitz near Graz, Austria.

14. Горбунов М. Е., Шмаков А. В., Среднестатистическая модель угла атмосферной рефракции на основе данных эксперимента COSMIC // Известия РАН, Физика атмосферы и океана, 2016, Т. 52, № 6, с. 699-706, doi: 10.7868/S0002351516060067.

15. Gorbunov, M. E., Ionospheric correction and statistical optimization of radio occultation data // Radio Sci., 2002, 37(5), 1084, doi:10.1029/2000RS002370.

16. Горбунов М.Е., Лауритсен К.Б., Родин А. Анализ данных эксперимента CHAMP по радиозатменному зондированию атмосферы Земли // Изв. РАН. Физика атмосферы и океана. 2005. V. 41. № 6, с. 798-813.

17. Gorbunov M.E., Kirchengast G. Wave-optics uncertainty propagation and regression-based bias model in GNSS radio occultation bending angle retrievals // Atmospheric Measurement Techniques, 2017, N 11, pp. 111-125. doi:10.5194/AMT-11-111-2018

18. Golub, G. H.; Reinsch, C. Singular value decomposition and least squares solutions // Linear Algebra. 1971, pp. 134-151. doi:10.1007/978-3-662-39778-7_10.

19. Golub, Gene H.; Van Loan, Charles F. An analysis of the total least squares problem // SIAM Journal on Numerical Analysis. 1980, 17 (6): 883-893. doi:10.1137/0717073.

20. Forsythe, G. E., Malcolm, M. A., Moler, C. B., Computer methods for mathematical computations. Englewood Cliffs, New Jersey 07632. Prentice Hall, Inc., 1977. XI, p.259