National software for processing information of the satellite – satellite space geodetic system

Full Text

1. ВВЕДЕНИЕ

17 марта 2002 г. произошло знаменательное событие, которое открыло новую эпоху в развитии космической геодезии: с российского космодрома Плесецк были выведены на околоземную близполярную орбиту два спутника-близнеца миссии GRACE с первоначальной высотой орбиты порядка 500 км и наклонением 89.0°, на расстоянии около 200 км друг от друга [1]. Близполярная орбита обеспечивала полное и относительно однородное покрытие земной поверхности измерениями. Начальная высота орбиты 500 км обеспечивала, с одной стороны, достаточный уровень гравитационного сигнала для уточнения ГПЗ, а с другой — подходила для исследований атмосферы и ионосферы. Спутники GRACE были связаны между собой высокоточной межспутниковой микроволновой системой, работающей в диапазонах К (24 ГГц) и Ка (32 ГГц) [2, 3]. Каждый космический аппарат (КА), помимо микроволновой системы измерений межспутникового расстояния, был оснащен GPS приемниками [4], звездными датчиками и прецизионными акселерометрами [5]. Высота орбиты спутников вследствие торможения в верхних слоях атмосферы уменьшалась в среднем на ≈ 30 м/сут, поэтому проекция трасс КА на поверхность Земли не была фиксированной и менялась со временем. В штатном режиме измерений положение спутников стабилизировалось по всем трем осям, поэтому антенны микроволновой системы удерживались точно направленными друг на друга. Микроволновый интерферометр К/Ka диапазона позволял измерять изменения межспутникового расстояния с микрометровой точностью. В сочетании с GPS измерениями, а также измерениями негравитационных ускорений и вращений, которые регистрировались бортовыми акселерометрами и звездными датчиками, микроволновая система, в конечном итоге, регистрировала воздействие вариаций ГПЗ на расстояние между спутниками.

Основной научной задачей миссии GRACE было создание ежемесячных моделей ГПЗ с пространственным разрешением порядка 400 км. За 15-летний срок активного существования миссия GRACE обеспечила мировое научное сообщество информацией о вариациях глобального гравитационного поля Земли и перемещениях больших масс во всех геосферах. Наблюдения GRACE легли в основу новых данных об изменениях климата [6], эволюции криосферы [7], океанов [8, 9, 10], гидрологических циклов [11, 12], а также тектонических событий [13].

Миссия GRACE доказала свою эффективность не только в получении моделей ГПЗ высокого пространственного разрешения, впервые по данным миссии было начато построение ежемесячных моделей. Поэтому в продолжение миссии GRACE 22 мая 2018 г. состоялся запуск двух спутников-близнецов GRACE Follow On. Основной целью миссии является обеспечение поступления данных для ежемесячных глобальных моделей ГПЗ высокого разрешения, начатых в миссии GRACE, т. е. мониторинг медленных временных вариаций ГПЗ. Параметры орбит спутников GRACE-FO практически те же самые, что и для спутников GRACE. С учетом опыта работы GRACE многие бортовые приборы и датчики были модернизированы [14]. Основным инструментом для измерения межспутникового расстояния по-прежнему является микроволновая система К-диапазона. Однако принципиальным моментом миссии является то, что спутники также оснащены лазерной интерферометрической системой, которая позволяет поднять точность измерений на нанометровый уровень. Первые данные обработки лазерных измерений превзошли все ожидания в отношении точности определения межспутникового расстояния [15], а также показали хорошее согласие с данными основной измерительной системы микроволнового диапазона [16].

С началом поступления данных GRACE, а с 2019 г. и данных GRACE-FO, в различных научных центрах мира начали создаваться модели ГПЗ (гравитационного поля Земли), как статические, так и ежемесячные. Полный список таких моделей ГПЗ, который регулярно пополняется, публикуется на сайте Международного центра глобальных моделей Земли (International Centre for Global Earth Models, ICGEM) [17]. Как правило, статические модели создаются на основе комбинации различных данных, которые могут включать данные космических геодезических миссий (CHAMP, GOCE, GRACE), данные альтиметрии и наземной гравиметрической съемки. В качестве примера можно указать одну из последних моделей ГПЗ — модель WHU-SWPU-GOGR2022S 2023 г., созданную китайскими учеными по данным GOCE и GRACE в разложении до степеней и порядков 300. Наиболее полной комбинированной глобальной моделью ГПЗ является модель XGM2019e, реализованная в разложениях 760, 2190 и 5540 сферических гармоник. Последнее соответствует пространственному разрешению около 4 км [18].

Тем не менее существуют статические модели ГПЗ, построенные только по данным миссии GRACE. Самой последней моделью, в которой были использованы все наблюдения GRACE, начиная с апреля 2002 г., является модель ITSG-Grace2018. Модель создана в Техническом университете города Грац (Австрия) и полна до степеней и порядков 200 [19, 20]. В 2016 г. были опубликованы две модели, созданные в университетах Китая Tongji (Шанхай) и HUST (Ухань): модель Tongji-Grace02s в разложении до 180 сферических гармоник [21] и модель HUST-Grace2016s в разложении до 160 сферических гармоник [22]. Глобальная модель ГПЗ ITU_GRACE16, основанная также только на данных миссии GRACE, разработана коллаборацией ряда научных организаций Турции в разложении до степеней и порядков 160 [23]. Можно также упомянуть более ранние модели ГПЗ, основанные только на данных GRACE. Это две модели, разработанные в Институте геодезии и геоинформации Университета Бонна: модель ITG-Grace2010s [24] и модель ITG-Grace03 [25]. Обе модели представлены в разложении до 180 сферических гармоник.

Очевидно, что модели ГПЗ, меняющиеся во времени на интервалах один месяц, могут быть построены только по данным GRACE и GRACE-FO. Такие модели создаются в трех международных центрах: Center for Space Research at University of Texas, Austin (CSR), Helmholtz Centre Potsdam German Research Centre for Geosciences (GFZ), Jet Propulsion Laboratory JPL). Решения, полученные в этих центрах, также доступны на сайте ICGEM [17].

В данной работе мы представляем первый отечественный программный комплекс, разработанный в ГАИШ МГУ, предназначенный для обработки данных наблюдений космической группировки типа GRACE, с целью решения обратной задачи по восстановлению гравитационного поля Земли и получения новой независимой модели ГПЗ.

2. ПРИНЦИП РАБОТЫ И ОСНОВНЫЕ КОМПОНЕНТЫ ПРОГРАММНОГО КОМПЛЕКСА

Программный комплекс обработки космической геодезической информации системы «спутник–спутник» предназначен для моделирования и решения целевых задач по восстановлению ГПЗ в среднечастотной области спектра. Термин «среднечастотный» означает, что рассматривается разложение ГПЗ в ряд по сферическим функциям до степени N_max ≤ 200 определенное выражением:

$R=\left(U-U_0\right)\approx \frac{GM}{r}\left[\begin{array}{l}\sum _{n=2}^{N_{\max }}\sum _{m=0}^n{\left(\frac{a_E}{r}\right)}^n{\overline{P}}_{nm}\left(\sin \phi \right)\times \\ \times \left({\overline{C}}_{nm}\cos m\lambda +{\overline{S}}_{nm}\sin m\lambda \right)\end{array}\right]$. (1)

В выражении (1): GM — гравитационный параметр Земли; r, φ, λ — сферические координаты точки в системе координат, жестко связанной с Землей (система ITRF); U₀ — потенциал Земли в представлении материальной точки, U₀ = GM/r ; a_E — средний экваториальный радиус Земли; ${\overline{P}}_{nm}\left(\sin \phi \right)$ — присоединенные функции Лежандра; ${\overline{C}}_{nm}$, ${\overline{S}}_{nm}$ — нормированные коэффициенты Стокса (m = 0 — зональные гармоники, m > 0 — тессеральные); n, m — степень и порядок разложения.

Блок-схема ПК, поясняющая его работу, представлена на рис. 1. В самом начале задаются параметры работы ПК (дата, количество восстанавливаемых гармоник ГПЗ, состав действующих сил на КА и т. д.), которые прописываются в отдельном текстовом файле. Этот файл считывается головной программой ПК, после чего идет обращение к файлам моделей, эфемерид, орбит КА (при наличии) и т. д., для выборки и записи в оперативную память ЭВМ необходимых данных для расчета. Условно эту часть ПК можно назвать модулем инициализации необходимых переменных.

 

Рис. 1. Блок-схема работы ПК

 

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

С учетом заданного состава сил формируется вектор правых частей дифференциальных уравнений движения КА, и посредством модуля интегратора рассчитываются необходимые величины системы «спутник–спутник», включая межспутниковые величины: координаты и компоненты скорости двух КА, межспутниковые расстояния, скорости и ускорения.

Рассчитанные величины позволяют составить матрицу условных и нормальных уравнений для дальнейшего восстановления ГПЗ методом наименьших квадратов (МНК).

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

Все воздействия сил на КА, а также алгоритм перехода из земной в небесную систему координат рассчитываются в ПК в соответствии со стандартами IERS2010 [26]. Состав используемых моделей приведен в табл. 1.

 

Таблица 1. Состав моделей сил, используемых в ПК

№	Тип модели	Наименования модели
1	Глобальная модель ГПЗ в виде сферических гармоник геопотенциала	EGM96 [27], ITG-Grace03, EGM2008 и другие с ресурса ICGEM [17]
2	Модель океанических приливов	FES2004/2014 [28, 29]
3	Модель приливов в теле Земли, включая полюсные приливы	Согласно IERS2010 [26]
4	Модель атмосферы	NRLMSISE-00 [30]
5	Модель гравитационного воздействия со стороны Луны, Солнца и планет	Согласно [31] на основе эфемерид DE421 [32]
6	Модель радиационного давления Солнца	Согласно [31], форма КА в виде прямоугольного параллелепипеда
7	Эффекты общей теории относительности	Согласно IERS2010 [26]

 

ПК имеет два основных режима работы:

• работа с синтетическими (модельными) данными;
• работа с реальными данными.

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

Далее запускается восстановление коэффициентов ГПЗ на основе рассчитанного файла. При восстановлении делается предположение, что все измерения можно распределить на группы так, чтобы внутри каждой группы измерения можно было считать равноточными. Можно предположить, что группа таких измерений может быть выполнена в течение некоторого времени, когда условия измерений не изменяются. Таким образом, все измерения группировки КА разбиваются на небольшие интервалы (дуги) продолжительностью 30 минут. После этого рассчитываются дуги и матрицы коэффициентов для восстановления ГПЗ. Коэффициенты записываются в отдельные файлы в рабочем каталоге программы. Каждая дуга состоит из набора положений, скоростей, ускорений. Для каждой точки дуги производится расчет относительных ускорений, рассчитанное значение сравнивается с измеренным и формируется матрица условных уравнений для коэффициентов ГПЗ.

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

После расчетов запускается непосредственно восстановление коэффициентов ГПЗ. Восстановление коэффициентов производится с помощью МНК с использованием методов библиотеки LAPACK [33] — библиотеки с открытым исходным кодом, содержащей методы для решения основных задач линейной алгебры. После выполнения МНК записывается новый файл коэффициентов ГПЗ. Запись коэффициентов идет до степени N_max, заданной в файле инициализации, в котором также можно задать выдачу разницы высот геоида получившегося и базового (начального) поля, а также дисперсию степеннóй разности, о которой более подробно будет сказано далее.

При работе с реальными данными вместо расчетных орбит и межспутниковых расстояний используются величины из оригинальных файлов миссий GRACE и GRACE-FO уровня 1B. В остальном режим аналогичен режиму с синтетическими данными.

3. СТРУКТУРА РЕАЛЬНЫХ ДАННЫХ

В качестве реальных данных были использованы данные миссий GRACE и GRACE-FO. Эти данные обладают схожей структурой (детальное описание структуры данных содержится в [34, 35]). Основной массив информации разделен на уровни, названные 1А и 1В. Уровень 1А содержит сырые данные, полученные со спутников, предварительно обработанные для приведения к читаемому стандартному виду. При этом учитывается калибровка приборов, в результате чего данные из бинарного формата переводятся в приборные значения и осуществляется их привязка к бортовой шкале времени. На уровне 1В данные профильтрованы и обработаны до вида, пригодного для восстановления гравитационного поля тем или иным методом. В частности, исключаются выбросы в показаниях приборов и все данные привязываются к единой шкале времени GPS с одинаковым шагом. Каждый уровень состоит из наборов посуточных файлов, в которых сгруппированы данные схожего типа.

Для решения основной задачи о восстановлении ГПЗ, в ПК были использованы следующие данные уровня 1B (приводятся аббревиатуры, с которых начинаются имена файлов данного типа): ACT1B — данные акселерометров; GNI1B и GNV1B — навигационные данные; KBR1B — данные измерения межспутникового расстояния микроволновой системой K/Ka диапазона; SCA1B — данные звездных датчиков.

Главный массив наблюдательных данных содержится в файлах типа KBR1B. Дальномерные микроволновые измерения, проводимые на борту с дискретностью 0.1 с., пересчитываются в межспутниковые расстояния с дискретностью 5 с. Этот массив данных содержит межспутниковое расстояние с точностью порядка 1 мкм, его первую и вторую производные по времени, световую и релятивистскую поправку к этим трем величинам, а также поправку за сдвиг фазового центра микроволнового излучателя относительно центра масс КА. В реализованном нами методе восстановления поля по межспутниковым ускорениям в качестве реальных данных использовались ускорения из KBR1B (фактически для каждого момента времени это сумма межспутникового ускорения, релятивистской поправки к нему и ускорения излучателя относительно центра масс).

Данные о кинематических орбитах спутников содержатся в файлах типов GNI1B и GNV1B. Отличаются эти типы системами отсчета, которые используются для представления данных. В GNI1B даны положения и скорости обоих аппаратов в земной системе отсчета GCRF (Geocentric Celestial Reference Frame), связанной с центром масс Земли и не вращающейся относительно барицентрической системы ICRF (International Celestial Reference Frame). В GNV1B положения и скорости даются во вращающейся вместе с Землей системе ITRF (International Terrestrial Reference Frame). В обоих типах файлов шаг данных составляет 1 с. Описанные орбитальные данные обоих типов использовались как для модельных вычислений, так и для вычислений по реальным данным.

Файлы типа ACT1B содержат данные акселерометров. Основная часть данных представляет собой негравитационные ускорения каждого спутника по показаниям бортовых акселерометров, которые приводятся с шагом 1 с. Негравитационные ускорения возникают главным образом в результате сопротивления атмосферы и воздействия солнечного излучения, прямого и отраженного от Земли. Нами были реализованы два режима учета такого рода эффектов: непосредственным вычислением по моделям негравитационных сил, а также путем использования данных акселерометров из файлов ACT1B.

Стоит отметить, что негравитационные ускорения в ACT1B приведены в системах отсчета SRF (Science Reference Frame), индивидуальных для каждого спутника, и требуют корректного пересчета в основную (невращающуюся) систему отсчета. Каждая такая спутниковая система жестко связана со своим спутником и устроена следующим образом. Начало отсчета помещено в центр масс пробного тела акселерометра (эта точка всегда находится в пределах 0.1 мм от центра масс спутника). Ось X направлена к фазовому центру микроволнового излучателя (по этой причине угол между вектором скорости летящего вторым спутника и его осью X всегда близок к 0°, а для летящего первым спутника этот угол близок к 180°). Ось Z перпендикулярна к оси X и к главной плоскости с оборудованием спутника, с направлением в сторону локального надира, то есть к поверхности Земли. Ось Y дополняет оси X и Z до правой тройки.

Преобразование от SRF к GCRF осуществляется с помощью информации от бортовых звездных датчиков, затабулированной в файлах типа SCA1B. По своей сути там дается ориентация каждого спутника (то есть направления осей его SRF) относительно земной невращающейся GCRF. В качестве способа задания такой ориентации выбрано представление в виде кватернионов. Как и орбитальные данные, так и данные акселерометров задаются с дискретностью 1 с. Каждый кватернион представляет собой набор из четырех чисел, с помощью которых по стандартным формулам (см., напр., [35]) делаются преобразования векторов от одной системы координат к другой (здесь от SRF к GCRF). При использовании в программе данных от акселерометров (ACT1B), так осуществляются необходимые повороты векторов негравитационных ускорений. В другом режиме, при вычислении модели негравитационных сил, данные файлов SCA1B применяются для нахождения ориентации отдельных элементов корпуса спутника относительно земной системы, что необходимо для корректного учета сопротивления воздуха и светового давления.

4. МАТЕМАТИЧЕСКИЕ ОСНОВЫ АЛГОРИТМА ВОССТАНОВЛЕНИЯ ГРАВИТАЦИОННОГО ПОЛЯ ЗЕМЛИ

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

$\begin{array}{c}{\ddot{r}}_i^o-{\ddot{r}}_i^c=\sum _{n=2}^{N_{\max }}\sum _{m=0}^n{a}_{nm}^i\Delta p_{nm},\\ a_{nm}^i={\left(\frac{\partial {\ddot{r}}^c}{\partial p_{nm}}\right)}_i,\\ p_{nm}=C_{nm},S_{nm},\end{array}$ (2)

где ${\ddot{r}}_i^o$ — измеренное взаимное ускорение между двумя КА, а ${\ddot{r}}_i^c$ — рассчитанное на основе модели исходного гравитационного поля с помощью алгоритма Беликова и Тайбаторова, используемого в ПК [36]; C_nm, S_nm — уточняемые коэффициенты ГПЗ. Выражение для расчета ${\ddot{r}}_i^c$ имеет следующий вид:

$\begin{array}{c}{\ddot{r}}^c={\ddot{r}}_{12}{e}_{12}+\frac{1}{D}{\dot{r}}_{12}{\dot{r}}_{12}-\frac{{\dot{D}}^2}{D},\\ {\ddot{r}}_{12}={\ddot{r}}_2-{\ddot{r}}_1,\\ {\dot{r}}_{12}={\dot{r}}_2-{\dot{r}}_1,\end{array}$ (3)

где $r_i={\{X_i,Y_i,Z_i\}}_{i=\mathrm{1,2}}$ — радиус-вектор каждого КА в инерциальной системе координат GCRS; ${\dot{r}}_i={\{V_{x_i},V_{y_i},V_{z_i}\}}_{i=\mathrm{1,2}}$ — вектор скорости каждого КА; ${\ddot{r}}_i={\{a_{x_i},a_{y_i},a_{z_i}\}}_{i=\mathrm{1,2}}$ — вектор ускорения каждого КА; $D=\sqrt{{(X_1-X_2)}^2+{(Y_1-Y_2)}^2+{(Z_1-Z_2)}^2}$ — расстояние между КА; $\dot{D}$ — скорость изменения расстояния между КА.

Для производных $a_{nm}^i$ получим:

$\begin{array}{c}{a}_{nm}^i={\left(\frac{\partial {\ddot{r}}^c}{\partial p_{nm}}\right)}_i=\frac{(r_1-r_2)}{D}\frac{\partial ({\ddot{r}}_1-{\ddot{r}}_2)}{\partial p_{nm}};\\ p_{nm}=C_{nm},S_{nm}.\end{array}$ (4)

Производные вектора ускорения по коэффициентам Стокса вычисляются на основе алгоритма Беликова и Тайбаторова.

Веса наблюдений w_i, (i = 1, 2, ..., n_obs), где n_obs — количество наблюдений, вводятся на этапе составления условных уравнений

$\begin{array}{c}\Delta {\xi }_i\sqrt{w_i}=\sum _{k=1}^{n_{def}}\sqrt{w_i}{a}_k^{\left(i\right)}\Delta x_k, i=\mathrm{1,2,...,}{n}_{obs};\\ a_k^{\left(i\right)}=\frac{\partial {\ddot{r}}_i^c}{\partial x_k}=\frac{r_1-r_2}{D}\frac{\partial ({\ddot{r}}_1-{\ddot{r}}_2)}{\partial x_k}, x_k=C_{nm},S_{nm};\\ \Delta {\xi }_i={\ddot{r}}_i^o-{\ddot{r}}_i^c,\\ \Delta x={\left(\Delta C_{nm}, \dots , \Delta C_{n_{\max }{m}_{\max }}, \Delta S_{nm}, \dots , \Delta S_{n_{\max }{m}_{\max }}\right)}^T,\end{array}$ (5)

где n_def — количество восстанавливаемых параметров, которое зависит от максимальной восстанавливаемой степени n_max следующим образом:

$n_{def}=n_{\max }^2+2n_{\max }-3$.

Ковариационная матрица ошибок наблюдений будет иметь равные диагональные элементы, если ошибки левых частей условных уравнений будут примерно одинаковыми для всех i. Сделать так можно было бы, если бы погрешности величин Δξ_i были бы известными. Достичь указанного свойства ковариационной матрицы можно, если положить:

$\sqrt{w_i}=\frac{1}{{\sigma }_i}$, (6)

где σ_i есть ошибка измерения с номером i.

Проблема в том, что обычно мы не знаем ошибок σ_i, поэтому приходится принимать какую-нибудь подходящую гипотезу относительно ошибок измерений. Именно поэтому все измерения разбиваются на дуги, о которых шла речь выше, для которых могут быть составлены нормальные уравнения:

$N_i\Delta x=b_i$, (7)

где

$\begin{array}{c}{N}_i=\left(\begin{array}{cccc}{l}_{11}& l_{12}& \dots & l_{1n}\\ l_{21}& l_{22}& \dots & l_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ l_{n1}& l_{n2}& \dots & l_{mn}\end{array}\right),\\ b_i={\left(b_1, b_2, \dots b_n\right)}^T,\\ l_{kj}=\sum _{i=1}^m{w}_i{a}_k^{\left(i\right)}{a}_j^{\left(i\right)}, k,j=\mathrm{1,} \mathrm{2,} \dots n,\\ b_k=\sum _{i=1}^m{w}_i{a}_k^{\left(i\right)}\Delta {\xi }_i, k=\mathrm{1,} \mathrm{2,} \dots n.\end{array}$ (8)

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

$N\Delta x=b$, (9)

где

$N=\sum _{k=1}^l\frac{1}{{\sigma }_k^2}{N}_k, b=\sum _{k=1}^l\frac{1}{{\sigma }_k^2}{b}_k$. (10)

Это возможно сделать в силу предположения, что ковариационные матрицы каждой из групп наблюдений независимы. Величина σ_k — ошибка измерений группы наблюдений k. После чего может быть найдено решение

$\Delta x=N^{-1}b$. (11)

Для анализа полученных результатов и сравнения различных моделей ГПЗ принято использовать разность полученного решения и эталонной модели в терминах высот геоида. Для разности высот геоида в представлении сферических функций используется выражение

$\begin{array}{c}\Delta N(\theta ,\lambda )=R_e\sum _{n=N_{\min }}^{N_{\max }}\sum _{m=0}^n\left[\Delta C_{nm}\cos \left(m\lambda \right)+\right.\\ \left.+\Delta S_{nm}\sin \left(m\lambda \right)\right]{\overline{P}}_{nm}\left(\cos \theta \right).\end{array}$ (12)

В выражении (12) величины ΔC_nm и ΔS_nm могут означать разности между коэффициентами Стокса разных моделей, либо коэффициентами Стокса полученного решения и эталонной модели. R_e — средний радиус Земли (опорный эллипсоид), n и m — степень и порядок разложения геопотенциала, ${\overline{P}}_{nm}\left(\cos \theta \right)$ — присоединенные функции Лежандра. Выражение (12) очень часто используется для представления результатов уточнения гравитационного поля Земли в силу своей наглядности. Но иногда бывает удобным описать точность решения по степеням гармоник, в этом случае дисперсия степенной разности (degree difference variance — DDV), выраженная в высотах геоида, есть

$DDV\left(n\right)=R_e\sqrt{\sum _{m=0}^n\left(\Delta C_{nm}^2+\Delta S_{nm}^2\right)}$. (13)

Оба этих похода к представлению результатов используются в ПК, однако далее в статье все результаты будут представлены в терминах дисперсии степеннóй разности.

5. РЕЗУЛЬТАТЫ ПРИМЕНЕНИЯ ПК К ВОССТАНОВЛЕНИЮ ПАРАМЕТРОВ ГПЗ ПО ДАННЫМ МИССИЙ GRACE И GRACE-FO

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

В настоящее время единственными данными для определения гравитационного поля Земли методом межспутниковых измерений являются ряды данных, полученные с 2002 г. по настоящее время в экспериментах GRACE и GRACE-FO.

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

В рамках настоящей работы были выбраны следующие эталонные поля.

1) ITSG-Grace2018s.gfc [37]: получено путем обработки измерений GRACE с 2002/04 до 2017/10, содержит коэффициенты до n = m = 200.

2) ITG-Grace2010s [37]: получено путем обработки измерений GRACE с 2002/08 до 2009/08, содержит коэффициенты до n = m = 180.

3) ITG-Grace03 [37]: получено путем обработки измерений GRACE с 2002/08 до 2007/04, содержит коэффициенты до n = m = 180.

4) Ежемесячные поля, публикуемые в рамках GRACE/GRACE-FO лабораториями GFZ, CSR, JPL (файлы Level 2 формата GSM), эти поля содержат коэффициенты до n = m = 60 или n = m = 96.

Все указанные эталонные поля получены только по данным GRACE и GRACE-FO без привлечения других спутниковых данных или наземных измерений.

5.1. Восстановление ГПЗ по данным GRACE на разных временных интервалах

На приведенных ниже рисунках, помимо дисперсии степенной разности восстановленного и эталонного гравитационных полей, приводится также кривая «geo», соответствующая кумулятивному гравитационному сигналу на степень (в терминах высоты геоида) эталонного ГПЗ без учета второй зональной гармоники C₂₀.

На рис. 2 представлен результат восстановления ГПЗ по данным GRACE на разных временных интервалах. В качестве эталонного поля использована модель ITSG-Grace2018s.gfc. Полученное решение на всех временных интервалах выглядит вполне достаточным, в среднечастотной части спектра (n = 40–50) разность между восстановленным и эталонным полями не превосходит долей миллиметра, в худшем случае (n ≈ 100) достигает единиц сантиметров.

 

Рис. 2. Результаты восстановления ГПЗ по данным GRACE до n = 100 в сравнении с моделью ITSG-Grace2018s.gfc

 

5.2. Восстановление ГПЗ по данным GRACE и GRACE-FO на 30-дневных интервалах

Ежемесячные модели GSM, рассчитанные по данным GRACE, содержат коэффициенты сферических гармоник до n = m = 90, а рассчитанные по данным GRACE-FO — до n = m = 96. В рамках тестовых расчетов выполнено восстановление для каждого из 12 месяцев 2010 и 2021 гг., сравнение осуществлялось с соответствующими моделями GSM (лаборатория GFZ).

Помесячное восстановление по данным GRACE (2010 г., рис. 3) показывает, что качество результата восстановления ПК достаточно стабильно, хотя в данных GRACE в целом довольно много пропусков измерений. В ПК реализован алгоритм исключения какой-либо дуги, если она содержит пропуск в данных. Если задать длину дуги в 2 часа, то, например, в марте 2010 г. отбрасывается 80% дуг, в мае — 74%, в июне — 63%; результат восстановления при этом становится ошибочен из-за недостаточности данных. Для результатов, приведенных на рис. 3, взята длина дуги равная 5 минутам.

 

Рис. 3. Результаты восстановления ГПЗ до n = 90 на 30-дневных интервалах по данным GRACE за 2010 г. в сравнении с моделями GSM (GFZ)

 

На рис. 4 представлены помесячные модели восстановленного ГПЗ по данным GRACE-FO за 2021 г. Качество данных у GRACE-FO выше, чем у GRACE, количество пропусков в данных значительно меньше. Соответственно, восстановленные модели ГПЗ лучше согласуются с моделями GSM, чем по данным GRACE. Разность между восстановленными моделями и моделями GSM в среднечастотной части спектра (n = 30–50) не превосходит 0.1–0.3 мм. Качество восстановления также стабильно.

 

Рис. 4. Результаты восстановления ГПЗ до n = 96 на 30-дневных интервалах по данным GRACE-FO за 2021 г. в сравнении с моделями GSM (GFZ)

 

5.3. Восстановление ГПЗ по данным GRACE-FO на длительных интервалах

В рамках тестовых расчетов было выполнено восстановление ГПЗ на длительных интервалах времени. Были использованы данные GRACE начиная с января 2003 г., поскольку качество данных в 2002 г., полученных в самом начале работы GRACE, заметно ниже и требует специальных подходов при обработке. Для сравнения были выбраны глобальные модели ITG-Grace03 и ITG-Grace2010s, которые были построены по данным GRACE, доступным на момент их разработки. Поэтому восстановление ГПЗ по данным миссии GRACE было выполнено на интервале 4.3 года, что соответствует модели ITG-Grace03, и на интервале 7.6 лет, что соответствует модели ITG-Grace2010s. Самой последней моделью, в которой были использованы все данные GRACE (15 лет наблюдений), является модель ITSG-Grace2018s. Однако восстановление гравитационного поля на интервале 15 лет нами не проводилось. Восстановление ГПЗ на указанных временных интервалах в сравнении с эталонными моделями ITG-Grace03 и ITG-Grace2010s в разложении до тех же степеней (n= 180) показано на рис. 5. Для сравнения на рис. 5 также показаны дисперсии степенных разностей моделей ITG-Grace2010s и XGM2019e, а также ITSG-Grace2018s и XGM2019e. Напомним, что модель XGM2019e является на сегодняшний день наиболее полной комбинированной глобальной моделью ГПЗ [18].

 

Рис. 5. Результаты восстановления ГПЗ до n = 180 по данным GRACE на интервалах 4.3 года и 7.6 лет в сравнении с моделями ITG-Grace03 и ITG-Grace2010s. Показано также сравнение моделей ITG-Grace2010s и ITSG-Grace2018s с комбинированной моделью XGM2019e

 

Дисперсии степенных разностей восстановленного поля и эталонных моделей демонстрируют характерное поведение ошибок коэффициентов сферических гармоник: сначала они уменьшаются, достигая минимального значения в низкочастотной и начале среднечастотной частей спектра, затем они возрастают, достигая соответствующих значений гравитационного сигнала (кривая «geo») (см., напр., [38]). Минимальных значений высот геоида (≈ 0.1 мм) эти разности достигают для степеней n = 20–60. Формально разложение имеет смысл, пока величина ошибки сферической гармоники не станет равной значению самой сферической гармоники. В данном случае это происходит при n = 160. Таким же образом меняется поведение дисперсий степенных разностей эталонных моделей в сравнении комбинированной моделью XGM2019e (кривые Grace2010s vs XGM2019e и Grace2018s vs XGM2019e). В целом изменение дисперсий степенных разностей восстановленного поля и эталонных моделей соответствует кривой Grace2010s vs XGM2019e, что говорит об адекватном восстановлении ГПЗ на длительных временных интервалах с использованием разработанного ПК.

6. ЗАКЛЮЧЕНИЕ

Космическая гравиметрия играет важную роль в изучении глобальных процессов, происходящих в различных геосферах Земли. Эти исследования опираются на модели гравитационного поля Земли, как статические, так и меняющиеся во времени, создаваемые в первую очередь по данным миссий GRACE и GRACE-FO. Моделей ГПЗ в настоящее время насчитывается несколько десятков. Активно в создании моделей ГПЗ участвуют различные университеты. Среди них два университета Китая — Tongji University (модель Tongji-Grace02s) и Huazhong University of Science and Technology (модель HUST-Grace2016s). Самой последней моделью, в которой были использованы все наблюдения GRACE, начиная с апреля 2002 г., является модель ITSG-Grace2018, созданная в Техническом университете города Грац (Австрия).

В ГАИШ МГУ разработан программный комплекс, позволяющий по данным космических миссий, аналогичных GRACE и GRACE-FO, восстанавливать параметры ГПЗ, и, таким образом, также строить различные модели ГПЗ. Результаты восстановления ГПЗ по данным миссий GRACE и GRACE-FO, полученные с использованием ПК ГАИШ МГУ, показали, что в целом решения соответствуют эталонным моделям как на длительных временных интервалах, так и на интервалах длительностью 1 месяц. На интервалах 4.3 года и 7.6 лет дисперсии степенных разностей восстановленного поля и эталонных моделей ITG-Grace03 и ITG-Grace2010s достигают минимальных значений высот геоида (≈ 0.1 м) для степеней n = 20–60. На 30-дневных интервалах разность высот геоида между восстановленными моделями и эталонными моделями GSM в среднечастотной части спектра (n = 30–50) не превосходит 0.1–0.3 мм.

ФИНАНСИРОВАНИЕ

Работа выполнена при финансовой поддержке Российского научного фонда, Грант № 23-42-00055.

About the authors

A. S. Zhamkov

Lomonosov Moscow State University, Sternberg Astronomical Institute

Email: zhamkov@physics.msu.ru
Russian Federation, Moscow

S. V. Ayukov

Lomonosov Moscow State University, Sternberg Astronomical Institute

Email: s.ayukov@gmail.com
Russian Federation, Moscow

A. I. Filetkin

Lomonosov Moscow State University, Sternberg Astronomical Institute

Email: alex-filyotkin@mail.ru
Russian Federation, Moscow

V. K. Milyukov

Lomonosov Moscow State University, Sternberg Astronomical Institute

Email: vmilyukov@yandex.ru
Russian Federation, Moscow

I. Yu. Vlasov

Lomonosov Moscow State University, Sternberg Astronomical Institute

Email: i.vlasov78@mail.ru
Russian Federation, Moscow

V. N. Sementsov

Lomonosov Moscow State University, Sternberg Astronomical Institute

Email: valera@sai.msu.ru
Russian Federation, Moscow

I. V. Gusev

Lomonosov Moscow State University, Sternberg Astronomical Institute

Email: GIV-WR@yandex.ru
Russian Federation, Moscow

V. E. Zharov

Lomonosov Moscow State University

Author for correspondence.
Email: vladzh2007@yandex.ru

Faculty of Physics

Russian Federation, Moscow

References

Supplementary files