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

УДК 525.2

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

ВРАЩЕНИЯ ЗЕМЛИ

JI. В. Зотов

(.ГАИШ) E-mail: tempus@sai.msu.ru

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

Введение

Идея о том, что наша планета может вращаться неравномерно, а ось вращения менять свое положение возникла еще в XVIII веке. Эйлером в 1765 г. была построена теория вращения твердой Земли, из которой следовало, что ось вращения может совершать круговые движения с периодом около 304 сут. Однако обнаружить периодические колебания оси вращения Земли удалось лишь в конце XIX в. Астроном Кюстнер, определявший постоянную аберрации, открыл годовой цикл изменения широт, а коммерсант Чандлер из Кембриджа, обработав ряды широтных наблюдений, выявил цикл с периодом около 430 сут, получивший его имя. Ньюкомб показал, что период Чандлера является аналогом периода Эйлера для Земли, обладающей упругими недрами и подвижным океаном. Эти открытия привели к созданию в 1899 г. Международной службы широты (МСШ). Сеть обсерваторий МСШ проводила наблюдения за изменениями широт, из которых определялось положение оси вращения Земли на протяжении более чем 80 лет.

Положение оси вращения в теле Земли задается двумя прямоугольными координатами X и У. Они отечитываютея от небесного эфемеридного полюса (НЭП), ось X направлена вдоль Гринвичского меридиана, ось У — в направлении 90° западной долготы. Ось, проходящая через НЭП, совершает прецессионные и нутационные колебания в инерци-альном пространстве. В 1962 г. задача определения координат полюса была возложена на Международную службу движения полюса (МСДП).

Предположение о возможном замедлении скорости вращения планеты под действием приливного трения высказал И. Кант в 1754 г. Подтвердилось это из наблюдения за движением планет Солнечной системы и анализа древних астрономических наблюдений лишь более века спустя. И только изобретение атомных часов во второй половине XX в. позволило начать регулярный и точный мониторинг отклонений в скорости вращения планеты. Сличе-

ние атомной шкалы Всемирного координированного времени UTC с получаемой из астрономических наблюдений и связанной с вращением Земли шкалой Всемирного времени UT1 позволяет отслеживать эти отклонения, а также вычислять продолжительность суток. В 1962 г. к работе приступило Международное бюро времени (МБВ), и стали публиковаться бюллетени, содержащие разность UT1 — UTC и величину продолжительности суток. Временной ряд UT1 — UTC отражает помимо скорости вращения Земли все изменения, которые претерпевала шкала UTC. Вследствие замедления вращения Земли шкала UT1 отстает от UTC. Шкала UTC используется повсеместно и время от времени подводится для согласования с вращением Земли. До 1972 г. к UTC при необходимости каждый раз добавлялась 0.1 с, так чтобы рассогласование не превышало 0.1 с. Однако столь частые поправки оказались неудобны. В 1972 г. было принято решение о том, что вносить изменения в UTC нужно реже и добавлять по 1 с, для того чтобы рассогласование с UT1 не превысило 0.7 с. Последняя такая секунда была добавлена к последней секунде 1998 г. и с тех пор, по причине ускорения вращения Земли, не добавлялась. По новым соглашениям следующая будет добавлена лишь тогда, когда рассогласование приблизится к 0.9 с, что может случиться в 2006 г.

В 1982 г. была создана Международная служба вращения Земли (МСВЗ), на которую были возложены задачи, ранее возлагавшиеся на МСДП и МБВ. Именно МСВЗ осуществляет сейчас сбор и обработку наблюдений современных обсерваторий, публикует бюллетени с ПВЗ, осуществляет прогноз. Современные средства наблюдений — радиоинтерферометры со сверхдлинными базами, системы глобального позиционирования GPS и Глонасс, станции лазерной локации Луны и искусственных спутников, объединенные в международные сети, дают высокоточные данные о вращении Земли с временным разрешением, достигающим нескольких часов и даже минут. Погрешности определения координат полюса

составляют сотни микросекунд дуги, а погрешности величин 1Л1 — иТС — десятки микросекунд времени. Однако координаты полюса X и ¥ и ход вращения Земли иТ1 моделируются с точностью, все еще уступающей точности наблюдений.

Вращение планеты — уникальный процесс, отражающий множество геофизических и астрономических явлений. На него влияет большое число факторов, среди которых вариации приливного потенциала, обусловленного действием небесных тел, изменения момента импульса ветров, течений, годовой цикл возбуждения атмосферы, явление Эль-Ниньо, таяние ледников и др. [1]. Так или иначе все явления, приводящие к перераспределению масс оболочек Земли и момента импульса между ними, влияют на вращение Земли. Прогноз многих из этих факторов краткосрочен, а некоторых, как, например, землетрясений, пока вообще невозможен. Даже природа самого мощного во вращении Земли чандлеровеко-го колебания, амплитуда которого достигает 0.3 с дуги (около 10 м на поверхности Земли), остается туманной. Обнаружены вариации его амплитуды и периода. Возбуждающая его сила до сих пор точно не известна [2]. Это сильно затрудняет прогнозирование.

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

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

Исходные данные

В качестве исходных данных использовались бюллетени ЕОРС01 и ЕОРС04 МСВЗ. Бюллетень ЕОРС01 содержит координаты полюса с 1846 по 1889 г. с шагом в 0.1 года и с 1890 г. с шагом в 0.5 года. Данные по скорости вращения Земли имеются с 1962 г. На рис. 1 представлен график изменения координат полюса с 1900 г. На рис. 2 представлен временной ряд 1ГП — иТС. Графики зависимости погрешности от эпохи наблюдений иллюстрируют колоссальный прогресс в точности наблюдений. Бюллетень ЕОРС04 содержит параметры вращения Земли с 1962 г. по настоящее время с шагом в одни сутки.

Методика вычислений

Двухмесячный прогноз строился по шестилетнему базовому отрезку ряда. На первом этапе выполнялось моделирование тренда. Для координат полюса X и У использовалась линейная модель. Временной ряд 1Л1 — иТС предварительно спрямлялся, устранялись скачки, обусловленные введением дополни-

2000-

1960 -

0 8 1920 —

-0.2 X, сек. дуги 0

0.4

0.2

0 У, сек. дуги

| I I I ШИ| ! ! !!1И1| ! I 11III 1| ГТТТШ| I ! I Ш1!|

0.01 0.1 1 10 100 1000 миллисек. дуги

Рис. 1. Изменения координат полюса с 1900 г. (слева) и погрешности измерений (справа)

100

10

0.1

0.01

0.001

1960 1970 1980 1990 2000 2010 годы

1960 1970 1980 1990 2000 2010 годы

Рис. 2. Расхождение шкал времени UT1 — UTC (слева) и погрешности измерений (справа)

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

Второй этап — моделирование гармонических составляющих рядов, освобожденных от тренда. Перед тем, как приступить к их моделированию и прогнозу был выполнен спектральный анализ. Спектрограммы координат полюса X, Y и UT1 — UTC вычислялись с использованием преобразования Фурье по данным бюллетеня EOPCOl (рис. 3). Известно, что основными гармоническими составляющими в движении полюса являются годовая составляющая и чандле-ровское колебание со средним периодом 435 сут, однако период и амплитуда последнего непостоянны [2]. В составе ряда UT1 — UTC ярко выражены 18.6-летняя, годовая и полугодовая составляющие. Для более подробного исследования спектрального состава рядов координат полюса X и Y был выполнен вейвлет-анализ, в котором использовался вейвлет Морле с параметром а = 100 [3]. Трехмерные екалограммы на рис. 4 иллюстрируют эволюцию периодических составляющих во времени. Видно, как сильно уменьшились период и амплитуда

чандлеровского колебания в 1930-е гг. На базовом отрезке, по которому выполнялся прогноз координат X или Y, параметры годовой и чандлеровской гармоники подбирались нелинейным МНК. Таким же образом для UT1 — UTC подбирались параметры годовой и полугодовой гармоник.

Поведение оставшихся после извлечения полиномиального и гармонического трендов временных рядов рассматривалось как случайное, и на третьем этапе для их моделирования использовались статистические методы. Метод авторегрессионного (АР) моделирования основан на представлении отсчетов временного ряда моделью

х% —

N

£

j=i

(X q X j —

■Щ,

где aj — параметры модели, щ — поступающий на вход белый шум, N — порядок модели. Параметры aj вычислялись по алгоритму Берга. Порядок модели N был принят равным 50, данный порядок приводит к достаточно хорошим значениям остаточной дисперсии, информационного критерия Акаике и вполне удовлетворительным прогнозам [4, 5].

1/год 1/год

Рис. 3. Спектрограммы координат полюса X к Y (слева) и UT1 — UTC (справа). По оси абсцисс отложены частоты

(число колебаний за год)

1990198019701960195019401930192019101900л

0.6

0.8

1

1.2

1.4

1.6

1.8

1990198019701960195019401930192019101900

0.6

0.8

1.2

1.4

1.6

1.8

Рис. 4. Скалограммы координат полюса X (слева) и У (справа). По оси абсцисс отложены частоты (число колебаний

за год), по вертикальной — годы

Метод ереднеквадратичеекой коллокации (СКК) позволяет прогнозировать отсчеты временного ряда с привлечением статистической информации, заключенной в автоковариационной функции (АКФ). Смещенная оценка АКФ вычисляется по временному ряду по формуле

ЯххЫ) =

1

N

Ы-т-1 п=0

Если нам известна АКФ для числа точек Ж, превышающего число точек N1 базового отрезка, по которому строится прогноз, мы можем построить прогноз на N^N1 точек вперед. Для этого значения АКФ заносятся в симметрическую ковариационную матрицу С^хх, так что = Д(|г — - Прогноз выполняется по формуле

/ = ЯпЯй1^

где (5и — левая верхняя часть матрицы СЦХх размерности N1 х N1, а Qft — левая нижняя часть матрицы СЦхх размерности (Ж — N1) х N1. Вектор прогноза / имеет размерность N — N1, вектор базового сигнала I — N1 [6].

Для прогнозирования ПВЗ предпринята попытка использования нейронных сетей (НС). На начальном

этапе была построена модель «амеба», состоящая из одного нейрона с линией задержки на входе. Схема такой модели представлена на рис. 5. Отсчеты временного ряда р{г) поступают на вход линии задержки Б размерности N = 100, с линии задержки на к-м шаге вектор рй{к) поступает на нейрон, где взвешивается с весами ги(г) и смещается на величину Ъ. После прохождения результирующего сигнала п(к) через линейную передаточную функцию формирует-

Г-N

-У

Ч_/

N

Рис. 5. Структурная схема нейронной сети «амеба»

Сравнение средней точности прогнозов, полученных разными методами для разных интервалов времени

Горизонт Прогноз МСВЗ АР-прогноз НС-прогноз СКК-прогноз

сутки Х,У 0.001" ии - итс 0.001 с х,У 0.001" ии - итс 0.001 с Х,У 0.001" ии - итс 0.001 с Х,У 0.001" ии - итс 0.001 с

1 0.15 0.02 0.15 0.02 0.14 0.034 0.27 0.032

5 3 1.5 1.0 0.4 0.78 0.39 1.8 0.43

10 6 2.1 2.4 1.1 3.3 0.96 4.0 1.07

30 12 3.5 11.7 7.2 10.6 5.7 13.4 6.0

0.160.12-

Е 0.08ч

I 0.040-0.04 -

-0.44 -

0.530.52 -0.51 -0.500.49 -0.480.47-

\

53200

53220 MJD

53240

53260

53200

53220 MJD

53240

53260

-0.45 -

-0.46 -

-0.47 -

UT1-UTC СКК прогноз АР прогноз прогноз НС контрольные данные Прогноз МСВЗ

53200

53220 MJD

53240

53260

Рис. 6. Двухмесячный прогноз рядов координат полюса X (слева), Y (справа) и UT1 — UTC (внизу). По оси абсцисс

отложены модифицированные юлианские даты

ся выход а(к), представляющий собой в нашей модели прогноз будущего значения ряда. Таким образом, N предыдущих значений используются для прогноза будущего значения, которое в свою очередь может быть вновь использовано на входе. Весовые коэффициенты нейрона адаптируются для решения задачи с помощью тренировки. Сигнал неоднократно пропускается через сеть, и прогнозное значение сравнивается с контрольным. К модели «амеба» можно привести нейронную сеть, содержащую любое количество линейных слоев и один нейрон выходного слоя.

Для вычисления средних отклонений прогнозных значений от наблюдений было сделано 20 прогнозов для разных эпох в прошлом и выполнено их сравнение с реальными данными. Результаты представлены в таблице. Там же представлены средние погрешности прогнозов МСВЗ [7-9]. Можно видеть, что прогноз НС оказался наиболее адекватным из всех (рис. 6). На протяжении двух месяцев его точность для координат полюса оказалась лучше точности МСВЗ. АР и СКК прогнозы превышают по точности прогноз МСВЗ на интервалах до 15 сут для UT1 — UTC и до 30 сут для координат полюса X я Y.

Заключение

Проведенное исследование показывает, что прогнозирование ПВЗ регрессионными методами АР и СКК дает хорошие результаты, по точности сравнимые с прогнозами МСВЗ и даже лучше на первых неделях, прогноз нейронной сети превышает по точности прогнозы МСВЗ на всем интервале прогнозирования. Перспективными видятся подходы

к прогнозированию временных рядов, основанные на совместном использовании нейронных сетей и вейв-лет-анализа. Результаты прогноза ПВЗ доступны в Интернете [10].

Автор благодарит научного руководителя В. Е. Жарова за предложения по работе, В. Г. Сурдина за предложенное название для модели нейронной сети, О. В. Баринову за помощь в разработке программ. Работа выполнена при финансовой поддержке РФФИ (грант 02-05-39004).

Литература

1. Сидоренков Н.С. Физика нестабильностей вращения Земли. М„ 2002.

2. Yatskiv Y. // IAU Colloquium 178 Cagliari, Sardinia, Italy, 27-30 September 1999. ASP Conference Series. 2000. 208. P. 383.

Витязев B.B. Вейвлет-анализ временных рядов. СПб., 2001.

Марпл С.Л. Цифровой спектральный анализ и его приложения. М., 1990.

Лукаш Ю.П. Адаптивные методы краткосрочного прогнозирования временных рядов. М., 2003.

6. Губанов B.C. Обобщенный метод наименьших квадратов. СПб., 1997.

7. IERS Annual Report 2002. Frankfurt am Main, 2003. Kosek W., Kalarus M. // Artificial satellites. 2003. 38, N 2. P. 41.

Kosek W. // Journees 2002. Bucharest, 25-28 September 2002. P. 125.

10. http://lnfml.sai.msu.ru/~tempus/pvz/predication/index.litml.

Поступила в редакцию 13.09.04