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

Полный текст

Методы моделирования движения полюса

Задача моделирования движения земных полюсов (под которыми понимаются точки пересечения мгновенной оси вращения Земли с ее поверхностью) представляет значительный научный интерес в области теоретической и небесной механики, геофизики и астрономии. Исследование движения северного полюса имеет как фундаментальное значение для решения задач геофизики и небесной механики, так и практическую важность. Координаты земного полюса используются, прежде всего, для вычисления матрицы перехода между земной и небесной системами координат [1, 2].

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

Для решения указанной задачи применяются различные методы математического моделирования. Сложный процесс колебаний полюса характеризуется составляющими с существенно различающимися частотными и амплитудными характеристиками. Математическое описание колебаний земного полюса осуществляется как с помощью детерминированных, так и стохастических или комбинированных подходов [2–17].

Детерминированные модели движения земного полюса основываются, как правило, на теоретических исследованиях Л. Эйлера, определившего 305-суточный период свободной нутации для недеформируемой Земли, и модели С. Чандлера (1891 г.), выявившего по результатам многочисленных наблюдений вариации широт обсерваторий с двумя периодическими компонентами движения полюса — 365 и 430–440 суток [3, 4].

Расхождение между чандлеровским периодом и периодом прецессии Эйлера (305 суток для недеформируемой фигуры Земли) исследовалось в работах С. Ньюкома, А. Пуанкаре, Г. Джеффриса, А. Лява, У. Манка, Г. Макдональда, Ф. А. Слудского, М. С. Молоденского (см. обзоры в [3, 4]). В исторической традиции движение земного полюса с периодом 430–440 суток принято называть свободной нутацией деформируемой Земли, или чандлеровскими колебаниями полюса.

Для описания свободной нутации деформируемой Земли также применяются методы аналитической механики движения деформируемого твердого тела относительно центра масс в переменных Андуайе, или в переменных «действие – угол» [5, 6].

Обзор ранних исследований влияния аддитивных стохастических возмущений, зависящих от времени, на колебания земного полюса представлен в [7]. В работе [8] содержится детальный анализ построения комбинированных стохастических моделей движения земного полюса. В частности, в [9, 10] на основе детерминированной модели движения земного полюса и данных Международной службы вращения Земли¹(МСВЗ)  разработана комбинированная небесномеханическая стохастическая корреляционная модель движения Земли относительно центра масс.

Сложная динамика гидросферы, включающая мелко- и крупномасштабные движения, а также широкий спектр турбулентных флуктуаций, свидетельствует о нелинейном характере момента сил диссипации. В [11] исследовано влияние нелинейных флуктуационно-диссипативных моментов сил на устойчивость режимов чандлеровских колебаний.

Однако в ранее разработанных моделях как период, так и амплитуда обоих колебаний (чандлеровского и годичного) предполагаются постоянными, несмотря на их известную изменчивость в определенных пределах [12, 13]. Современная проблема переменности параметров основных компонент колебаний земного полюса остается актуальной и недостаточно изученной [14].

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

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

1. Математическая модель и перспективы ее развития

Как известно, простейшая модель движения земного полюса описывает двухчастотный колебательный процесс [16, 17]:
\[\begin{equation}
\tag{1}
x_{p} =c_{x} +a_{\rm ch} \cos w_{\rm ch} +a_{\rm h} \cos w_{\rm h} , 
\quad
y_{p} =c_{y} +a_{\rm ch} \sin w_{\rm ch} +a_{\rm h} \sin w_{\rm h} . 
\end{equation}\]
Здесь $c_{x}$, $c_{y}$ — координаты тренда, содержащие колебания с периодами свыше 6.45 лет; $a_{\rm ch}$, $a_{\rm h}$ и $w_{\rm ch}$, $w_{\rm h}$ — амплитуды и фазы чандлеровского и годичного колебаний соответственно.

При этом модель является двухчастотной в приближенном смысле — когда амплитуды и фазы основных составляющих постоянны. Применение такой двухчастотной модели имеет ограничения как по точности прогноза, так и по времени ее автономного использования. Как показано в [14, 18], на основании анализа длительных рядов наблюдений, амплитуды и фазы основных составляющих можно считать постоянными лишь на относительно небольших временных интервалах (до 20 лет). На таких интервалах параметры могут рассматриваться как квазипостоянные, определяемые на промежутках аппроксимации движения полюса, непосредственно предшествующих интервалам прогнозирования. В этом случае значения амплитуд и фаз основных составляющих колебаний полюса следует интерпретировать как средние величины, поскольку помимо медленных изменений на временных масштабах 40 лет и более эти параметры подвержены также незначительным вариациям с меньшими периодами.

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

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

2. Вариации параметров движения земного полюса с частотой прецессии орбиты Луны

В работе [14] предложено преобразование координат земного полюса, позволяющее выявить синфазность вариаций его движения с прецессией лунной орбиты. Переход от исходной земной системы координат $(x_{p}, y_{p})$ к новой системе $(\xi_{p}, \eta_{p})$, в которой полюс совершает колебания, синфазные с прецессионным движением лунной орбиты, в матричной форме описывается выражением
\[\begin{gather}
\tag{2}
\begin{pmatrix} 
{\xi_{p}} \\ {\eta_{p}} 
\end{pmatrix} 
= 
\Pi(w_{2}-w_{1}) 
\bigg[\Pi(w_{1}) 
\begin{pmatrix} {x_{p}-c_{x}} \\ {y_{p}-c_{y}} 
\end{pmatrix} - \begin{pmatrix} {a_{0}} \\ {0} \end{pmatrix}\bigg], 
\\ 
w_{2} = \begin{cases} 
w_{\rm h}, & \text{если  } a_{\rm h} < a_{\rm ch}, \\ 
w_{\rm ch}, & \text{если  } a_{\rm ch} < a_{\rm h}; 
\end{cases} 
\quad  
w_{1} = \begin{cases} 
w_{\rm ch}, & \text{если  } a_{\rm h} < a_{\rm ch}, \\ 
w_{\rm h}, & \text{если  } a_{\rm ch} < a_{\rm h}; 
\end{cases} 
\\
\dot{w}_{\rm h} = \nu\omega_{*}, \quad \dot{w}_{\rm ch} = N\omega_{*}. 
\end{gather}\]
Здесь $\Pi(\alpha)$ — матрица плоского поворота на угол $\alpha$; $a_{0}$ — среднее значение амплитуды колебаний полюса при его движении вокруг «средней точки» за 6-летний цикл (без учета трендовой составляющей); $c_{x}$, $c_{y}$ определяют положение «средней точки» полюса и содержат константы, вековые слагаемые и вариации с периодами более шести лет; $a_{\rm ch}$, $a_{\rm h}$ — амплитуды чандлеровской и годичной гармоник с фазами $w_{\rm ch}$, $w_{\rm h}$ соответственно; $N\cong 0.843$, $\nu = 1$ — чандлеровская и годичная частоты, измеряемые в циклах/год; $\omega_{*}$ — среднее движение барицентра системы «Земля – Луна» по орбите вокруг Солнца; $\dot{w}_{2} - \dot{w}_{1} = \pm \nu_{T}\omega_{*}$ — частота шестилетней цикличности движения полюса.

Отметим, что преобразование (2) представляет собой композицию преобразований «поворот» и «сдвиг», определяемую исключительно средними параметрами чандлеровской и годичной компонент. Как показано в [14], в системе координат $(\xi_{p},\eta_{p})$ полюс совершает колебания, синфазные с колебаниями угла наклона плоскости лунной орбиты к земному экватору. При этом колебания координаты $\xi_{p}$ синфазны с колебаниями угла наклона плоскости орбиты Луны к земному экватору, а фаза колебаний координаты $\eta_{p}$ сдвинута на $\pi/2$.

В настоящей работе ставится задача учета выявленных ранее вариаций параметров колебательного движения земного полюса с периодом 18.6 года (соответствующим периоду прецессии лунной орбиты и устанавливаемых посредством преобразования (2)) в модели вида (1). Кроме того, требуется оценить точность полученной расширенной модели с учетом дополнительных слагаемых, обусловленных долгопериодическим влиянием Луны.

3. Вариации амплитуд и фаз чандлеровской и годичной компонент в колебательном процессе земного полюса

Примем следующие допущения: положим $\dot{a}_{\rm ch} = 0$, $\dot{a}_{\rm h} = 0$, $w_{\rm ch} = Nt + \alpha_{\rm ch}$, $w_{\rm h} = \nu t + \alpha_{\rm h}$, где $\alpha_{\rm ch}$, $\alpha_{\rm h}$ — постоянные фазы. Таким образом, амплитуды $a_{\rm ch}$, $a_{\rm h}$ будем считать квазипостоянными, а фазы $w_{\rm ch}$, $w_{\rm h}$ — линейными функциями времени. Данные упрощения означают, что колебательный процесс земного полюса разделяется на тренд, чандлеровское и годичное колебания с квазипостоянными амплитудами и линейными фазами, при этом менее значимые гармоники исключаются из рассмотрения.

Подставляя (1) в (2), преобразуем компоненты вектора
\[\begin{equation}
\begin{pmatrix} 
{x_{1}} \\ {y_{1}} 
\end{pmatrix} = 
\Pi(w_{1})
\begin{pmatrix} {x_{p}-c_{x}} 
\\ {y_{p}-c_{y}} 
\end{pmatrix}
\end{equation}\]
к следующему виду:
\[\begin{multline*}
x_{1} = \frac{a_{\rm h} + a_{\rm ch}}{2} + 
\operatorname{sgn}\Bigl(\frac{a_{\rm h}}{a_{\rm ch}}\Bigr)\frac{a_{h} - a_{ch}}{2} + {}
\\
{}+ \Bigl[\frac{a_{\rm h} + a_{\rm ch}}{2} - \operatorname{sgn}\Bigl(\frac{a_{\rm h}}{a_{\rm ch}}\Bigr)\frac{a_{\rm h} - a_{\rm ch}}{2}\Bigr]\cos(w_{\rm h} - w_{\rm ch}), 
\end{multline*}\]
\[
y_{1} = \Bigl[\frac{a_{\rm h} - a_{\rm ch}}{2} - 
\operatorname{sgn}\Bigl(\frac{a_{\rm h}}{a_{\rm ch}}\Bigr)\frac{a_{\rm h} + a_{\rm ch}}{2}\Bigr]\sin(w_{\rm h} - w_{\rm ch}). \hspace{4cm}
\]

Приняв
\[
a_{0} = \frac{a_{\rm h} + a_{\rm ch}}{2} + 
\operatorname{sgn}\Bigl(\frac{a_{h}}{a_{\rm ch}}\Bigr)\frac{a_{\rm h} - a_{\rm ch}}{2},
\]
из (2) получаем
\[\begin{equation}
\tag{3}
\begin{pmatrix}{\xi_{p}} 
\\ {\eta_{p}} 
\end{pmatrix} = \begin{pmatrix} {\frac{a_{\rm h} + a_{\rm ch}}{2} - 
\operatorname{sgn}\bigl(\frac{a_{\rm h}}{a_{\rm ch}}\bigr)\frac{a_{\rm h} - a_{\rm ch}}{2}} 
\\ {0} 
\end{pmatrix}.
\end{equation}\]
Из (3) следует, что $\xi_{p} = a_{\rm h}$ при $a_{\rm h} < a_{\rm ch}$ и $\xi_{p} = a_{\rm ch}$ при $a_{\rm h} > a_{\rm ch}$.

Отметим, что для построения ряда $\xi_{p}$ по данным $x_{p}$, $y_{p}$ не требуется разделения колебаний на чандлеровскую и годичную компоненты, что также продемонстрировано в работе [14].

На основании приведенных выкладок можно построить преобразование вида (2) следующим образом:
\[\begin{equation}
\tag{4}
\begin{pmatrix} 
{\xi'_{p}} \\ {\eta'_{p}} 
\end{pmatrix} 
= \Pi(w_{1} - w_{2})
\bigg[
\Pi(w_{2})
\begin{pmatrix} 
{x_{p} - c_{x}} \\ {y_{p} - c_{y}} 
\end{pmatrix} - 
\begin{pmatrix} 
{a_{0}} \\ {0} 
\end{pmatrix}\bigg],
\end{equation}\]
при этом $\xi'_{p} = a_{\rm ch}$ при $a_{\rm h} < a_{\rm ch}$ и $\xi'_{p} = a_{\rm h}$ при $a_{\rm h} > a_{\rm ch}$.

Используя данные ряда C01 МСВЗ² о движении земного полюса за период 1900–2020 гг., построим ряд $\xi_{p} + \xi'_{p}$ с помощью (2) и (4). На рис.  1  представлен амплитудный спектр этого ряда. Значительная продолжительность интервала наблюдений позволяет выделить низкочастотные гармоники в окрестности рассматриваемой частоты. Наличие в спектре гармоники с периодом 18.6 года свидетельствует о том, что по крайней мере одна из амплитуд $a_{\rm ch}$, $a_{\rm h}$ чандлеровской и годичной компонент содержит гармонику с периодом 18.6 года. Для точного установления этого факта необходимо выделить чандлеровское и годичное колебания в движении полюса.

 

Рис. 1. Амплитудный спектр ряда \(\xi _{p} +\xi '_{p}\)
[Figure 1. The amplitude spectrum of the series \(\xi _{p} +\xi '_{p}\)]

 

Разделим спектр колебаний земного полюса на области. Для чандлеровской компоненты по каждой из координат выделим окрестность частоты 0.843 циклов/год с границами $0.843 \pm 0.157/2$ циклов/год. Аналогично, для годичной компоненты определим окрестность частоты 1 цикл/год с границами $1 \pm 0.157/2$ циклов/год. Совокупность чандлеровской и годичной компонент будет содержать все колебания из спектрального интервала $0.765\div1.0785$ циклов/год. Хотя такое разделение носит формальный характер, оно оправдано тем, что дополнительные гармоники, вызывающие 18-летнюю модуляцию основных компонент с частотой прецессии лунной орбиты, попадают в соответствующие области.

Для каждой компоненты можно построить ряды их амплитуд $a_{\rm ch}$, $a_{\rm h}$ и вариаций фаз $\delta w_{\rm ch}$, $\delta w_{\rm h}$ (в отличие от рассмотренного выше упрощенного случая линейных фаз, наблюдаемые фазы подвержены изменениям). Применение различных методов выделения гармоник приводит к согласующимся оценкам фаз и относительно небольшому разбросу в оценке средних амплитуд. Это позволяет аппроксимировать рассматриваемые 18-летние колебания параметров чандлеровской и годичной компонент стационарными гармониками, которые можно интерпретировать как гармоники со средними амплитудами.

На рис. 2 и 3 представлены стационарные гармоники со средними амплитудами, выделенные из спектров указанных параметров. Вариации амплитуд $a_{\rm ch}$, $a_{\rm h}$ чандлеровской и годичной компонент синфазны, но различаются по величине, тогда как вариации полярных углов $w_{\rm ch}$, $w_{\rm h}$ совпадают как по фазе, так и по амплитуде.

 

Рис. 2. Вариации амплитуд чандлеровской и годичной компонент с периодом 18.6 года прецессии лунной орбиты (в угловых миллисекундах); по оси абсцисс — время (годы)
[Figure 2. Variations in the amplitudes of the Chandler  and annual wobbles with the 18.6-year lunar orbital precession period (in milliarcseconds); the x-axis is time (years)]

 

Рис. 3. Вариации фаз чандлеровской и годичной компонент с периодом 18.6 года прецессии лунной орбиты (в радианах); по оси абсцисс — время (годы)
[Figure 3. Variations in the phases of the Chandler and annual wobbles with the 18.6-year lunar orbital precession period  (in radians); the x-axis is time (years)]

 

4. Численное моделирование колебаний земного полюса с учетом вариаций параметров чандлеровской и годичной компонент

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

Для этого необходимо аппроксимировать 18-летнюю цикличность параметров $a_{\rm ch}$, $w_{\rm ch}$, $a_{\rm h}$, $w_{\rm h}$ гармониками $\cos \Omega$ и $-\sin \Omega$, где выражение для $\Omega$ имеет вид [1, p. 67]
\[
\Omega \approx 125.04455501^{\circ} - 6962890.5431''t + 7.4722''t^{2} .
\]
Здесь $t$ — время в столетиях, отсчитываемое от 12 ч. 1 января 2000 года. После определения амплитуд этих гармоник выполним обратное преобразование для получения дополнительных слагаемых к модели (1) движения земного полюса:
\[\begin{gather}
\begin{pmatrix} 
{\Delta x_{\rm ch}} \\ {\Delta y_{\rm ch}} \end{pmatrix} = 
\Pi^{-1}(w_{\rm ch})
\begin{pmatrix} {a_{\rm ch} \cos \delta w_{\rm ch} - a_{\rm ch}^{0}} \\ {a_{\rm ch} \sin \delta w_{\rm ch}} \end{pmatrix}, \\
\begin{pmatrix} {\Delta x_{\rm h}} \\ {\Delta y_{\rm h}} \end{pmatrix} = 
\Pi^{-1}(w_{\rm h})\begin{pmatrix} {a_{\rm h} \cos \delta w_{\rm h} - a_{\rm h}^{0}} \\ 
{a_{\rm h} \sin \delta w_{\rm h}} \end{pmatrix},
\end{gather}\]
где $a_{\rm ch}^{0}$, $a_{\rm h}^{0}$ — средние значения амплитуд чандлеровской и годичной компонент соответственно. Их оценки могут быть получены, например, методом наименьших квадратов при аппроксимации (1) на длительном интервале наблюдений (с 1900 по 2020 гг.).

Модель колебаний полюса с учетом дополнительных слагаемых представляется в виде
\[\begin{equation}
\tag{5}
\tilde{x}_{p} = x_{p} + \Delta x_{\rm ch} + \Delta x_{\rm h}, \quad 
\tilde{y}_{p} = y_{p} + \Delta y_{\rm ch} + \Delta y_{\rm h}.
\end{equation}\]

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

Сравним точность экстраполяций моделей (1) и (5) с дополнительными слагаемыми. Верификация модели проводилась в автономном режиме без коррекции параметров.

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

Была проведена серия расчетов со сдвигом тестового интервала по временной шкале. Для определения параметров двухчастотной модели использовался фиксированный тестовый интервал длительностью 15 лет. Расчеты выполнялись итерационно со сдвигом начала тестового интервала с 1900 года до 2005 года с шагом 2 года. Как показали результаты, изменение длительности тестового интервала и увеличение шага не приводят к качественным различиям в результатах.

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

В численном эксперименте рассчитывались СКО двухлетних экстраполяций. Двухлетний интервал экстраполяции сдвигался с шагом 1 год от 1962 года до 2018 года. Каждому тестовому интервалу соответствовало среднее значение СКО двухлетних экстраполяций за указанный период. На рис. 4 представлена разность $\delta_{\sigma}$ средних СКО экстраполяций для моделей (1) и (5). Каждое значение среднего СКО отнесено к началу соответствующего тестового интервала.

 

Рис. 4. Разность $\delta_{\sigma}$ средних квадратических отклонений экстраполяций для моделей (1) и (5) (в угловых миллисекундах); по оси абсцисс — время (годы)
[Figure 4. The difference $\delta_{\sigma}$ between the root-mean-square errors of extrapolations using models (1) and (5) (in milliarcseconds); the x-axis is time (years)]

 

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

Заключение

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

Конкурирующие интересы. У нас нет конфликта интересов в отношении авторства и публикации этой статьи.
Авторский вклад и ответственность. Все авторы принимали участие в разработке концепции статьи и в написании рукописи. Авторы несут полную ответственность за предоставление окончательной рукописи в печать. Окончательная версия рукописи была одобрена всеми авторами.
Финансирование. Исследование выполнялось без финансирования.

---

¹https://www.iers.org/
²https://datacenter.iers.org/products/eop/long-term/c01/

Об авторах

Сергей Сергеевич Крылов

Московский авиационный институт (национальный исследовательский университет)

Автор, ответственный за переписку.
Email: compgra@yandex.ru
ORCID iD: 0000-0003-3267-6411
SPIN-код: 8634-7203
Scopus Author ID: 55453093500
https://www.mathnet.ru/rus/person232488

кандидат физико-математических наук, доцент; заведующий кафедрой; каф. 806 вычислительной математики и программирования

Россия, 125993, Москва, Волоколамское шоссе, 4

Зо Аунг Мьо

Московский авиационный институт (национальный исследовательский университет)

Email: myozawaung53@gmail.com
ORCID iD: 0009-0009-8557-5965
https://www.mathnet.ru/rus/person232487

аспирант

Россия, 125993, Москва, Волоколамское шоссе, 4

Вадим Владимирович Перепелкин

Московский авиационный институт (национальный исследовательский университет)

Email: vadim802@gmail.com
ORCID iD: 0000-0002-9061-4991
Scopus Author ID: 8263058800
ResearcherId: S-6900-2019
https://www.mathnet.ru/rus/person68736

доктор физико-математических наук; профессор; каф. 802 мехатроники и теоретической механики

Россия, 125993, Москва, Волоколамское шоссе, 4

Список литературы

Дополнительные файлы