CC BY
129
Геодезия
УДК 528.34:629.783 Ю.В. Сурнин СГГ А, Новосибирск
АДАПТИВНЫЙ АЛГОРИТМ АППРОКСИМАЦИИ МОДЕЛИ ВРАЩЕНИЯ ЗЕМЛИ
Рассматривается адаптивная аппроксимация в соответствии с требуемой точностью классической модели вращения Земли. Адаптивная модель настраивается на заданную пользователем точность с помощью замены двенадцати углов вращения тремя углами Кардано и чебышевского альтернанса. Значительное уменьшение числа тригонометрических членов повышает быстродействие компьютерных расчетов и позволяет использовать модель для решения различных задач космической геодезии, связанных с численным интегрированием дифференциальных уравнений движения спутников.
модель вращения Земли, адаптация, углы Кардано, чебышевский альтернанс.
Yu.V. Surnin SSGA, Novosibirsk
ADAPTIVE ALGORITHM
OF EARTH ROTATION MODEL APPROXIMATION
Adaptive approximation according to required accuracy of classical model of Earth rotation is considered. The adaptive model is adjusted with prescribed by user accuracy with the help of replacing of twelve rotation angles by three angles of Cardano and Chebyshev's alternance. Considerable diminution of trigonometric terms number increases the speed of computation and allows to use the model for solving different tasks of space geodesy, connected with numerical integration of differential equation of sputnic's motion.
Earth rotation model, adaptation, Cardano angles, Chebyshev alternance
Математическая модель вращения Земли в рамках Ньютоновой механики обычно представляется в виде произведения элементарных матриц поворота [1, 2, 4], которые учитывают прецессионно-нутационное движение оси мира относительно инерциального пространства, сидерическое вращение Земли, полярное движение мгновенной оси в теле Земли и неравномерность ее вращения. Чтобы уменьшить влияние погрешностей модели вращения Земли на взаимное преобразование земных и небесных координат, целесообразно в качестве инерциальной системы координат (ИСК) - точнее квазиинерциальной -выбрать небесную систему отсчета [5, 6], определяемую истинным равноденствием и истинным экватором фиксированной эпохи Tm, лежащей примерно в середине некоторого интервала [Tn, Tk]. Например, Tn, Tk - эпохи начала и конца какой-либо наблюдательной кампании.
Тогда матрица W направляющих косинусов осей общеземной системы координат (ОЗСК) относительно осей ИСК может быть описана двенадцатью элементарными поворотами:
63
Геодезия
R1 ( s0m Asm) R3( Ащm) R1 (s0m) R3(ZAm) R2( @Am) R3(zAm) R1 ( s0) X
x R3(Ащ) Ri(so +As) R3(-S) R2(xp) Ri(yp). (1)
В формуле Aym, Ащ, и Asm, As - компоненты астрономической нутации в долготе и наклоне соответственно в эпохи Tm и T; ZAm, 0Am, zAm - экваториальные параметры прецессии на интервале времени [Tm, T]; s0, и s0m - средние наклоны эклиптики к экватору в эпохи T и Tm; S - истинное звездное время на меридиане Гринвича на момент времени T; xp, yp - параметры полярного движения оси вращения в теле Земли в эпоху Т. Параметр неравномерности вращения Земли A UT1(T) включен в угол S.
Если вместо ньтоновой механики для модели вращения Земли использовать общую теорию относительности [1, 4], то дальнейшие преобразования матрицы W принципиально не изменятся. Изменится только сама исходная матрица W в выражении
У = Wx, (2)
которое осуществляет связь прямоугольных координат произвольного вектора x = {x1, x2, x3}, заданного в ОЗСК, с его координатами y = {y1, y2, y3} в ИСК.
Современная аналитическая модель вращения Земли [1, 2, 4] содержит более двух тысяч полиномиальных и тригонометрических членов. Синусы и косинусы значительно увеличивают время численного интегрирования уравнений движения искусственных спутников Земли в динамическом методе космической геодезии [6]. С целью увеличения скорости расчетов предлагается преобразовать традиционную формулу (1) к виду
W = R3(-O)0. (3)
В правой части формулы (3) выделен главный член R3(-0, содержащий некоторый аналог гринвичского звездного угла в [1, 4]. Угол в вводится, как функция земного времени TT (практически реализуемого атомными стандартами частоты), по формулам:
в = dm + rn(T - Tm); dm = во + rn(Tm - To); ю = 2п(1 + ^) радиан/сутки; (4)
во = 2п • 0,7790572732640 радиан; ц = 0,002737811911354480,
где ю - номинальная угловая скорость вращения Земли; To = 2451545,0 - стандартная эпоха J2000,0. Заметим, что в [4] для вычисления звездного угла в рекомендуется всемирное время UT1, здесь же используется земное время TT.
Остальная часть в формуле (3) - матрица 0 - описывает небольшие колебательные движения оси вращения Земли на интервале времени [Tn, Tk] выражением
0 = R3(6)W.
(5)
64
Геодезия
Матрицу вращения 0 можно представить двояко. C одной стороны, как матрицу направляющих косинусов
0 = [Г]' к'],
(6)
строками которой являются орты i', j', к' осей ИСК в проекциях на орты ОЗСК
i = {1, 0, 0}, j = {0, 1, 0}, к = {0, 0, 1}.
С другой стороны, эту же матрицу 0 можно представить в виде произведения трех матриц элементарных поворотов либо на углы Эйлера, либо на углы Кардано. При малых колебаниях экватора (по наклону) возникает квазилинейная зависимость двух углов Эйлера (прецессии и чистого вращения), примыкающих друг к другу по линии узлов. В связи с этим, отдавая предпочтение углам Кардано, правую часть формулы (6) представим как
0=R3(-A0)R2(e)R1 (о), (7)
где кардановые углы можно интерпретировать следующим образом: АО - аналог нутации по прямому восхождению; в - аналог нутации по склонению и о -аналог нутации по наклону. Углы Кардано АО, в, о можно вычислить через исходные орты i', j', к'. Последние получаются, как столбцы матрицы 0, по формулам (1), (5) и (6).
Если ввести три вспомогательных орта i", j", к" равенствами
j" = к'Щ к'М |; к" = ixj"; i" = ]"хк\ (8)
то в области определения [-п/2, п/2] углы Кардано найдутся по формулам:
tgA6 = -i'j"/i4"; tgfi = кЧ/к'к"; tgo = ку^к-к". (9)
Следует обратить внимание на то, что эмпирические углы xp, yp и (1 + n)AUT1 даются не аналитическими выражениями, как это имеет место для всех остальных углов прецессии и нутации в равенстве (1), а публикуются в виде таблицы по аргументу Т. Поэтому табличное задание трех функций xp(T), yp(T), (1 + /u)AUTl(T) лучше преобразовать в аналитический вид. Для этой цели, чтобы не снизить табличную точность, можно использовать, например, аппроксимацию кубическими сплайнами со сглаживанием [7].
Теперь малые углы Кардано A6(t), в(0, o(t) можно рассматривать как аналитические функции времени t, даваемые формулами (1), (5), (6), (8), (9) и аппроксимировать их на интервале [Ти, Tk] полиномами Чебышева, используя че-бышевский альтернанс и равномерную сходимость ряда [8].
Пусть f(x) любой из углов Кардано A6(t), e(t), o(t), t = T- Tm. Тогда разложение функции f(x) по полиномам Чебышева Tk(x) (k = 0, 1, ..., N - 1) до степени N - 1 можно записать так:
65
Геодезия
N-1
f (x) = Z ckTk(xX (10)
к=0 где
x = (t - to)/At; to = (tk + tn)/2; At = (tk - tn)2; t = T- Tm; -1 < x < 1;
To(x) = 1; Tj (x) = x; Tk+i (x) = 2xTk(x) - Ты(х); к = 1, N - 1.
Коэффициенты ck разложения (10) находятся по формулам:
j=N
Ck= [(2 - dok)/N\ Z f(xj)Tk(xj), k = 0, 1, 2, ..., N - 1, (11)
j=1
S0k = 1, если k = 0 иначе 0,
где функции f(xj) и полиномы Tk(xj) вычисляются с помощью равенств (1), (5), (6), (8), (9) в корнях xj полинома Чебышева TN степени N. Заметим, что cN = 0, поэтому в ряде (10) предел суммирования равен N - 1.
Для всех трех функций A6(t), fi(t), a(t) узлы xj вычисляются один раз (с минимальным обращением к функциям cos и sin) по рекуррентным формулам:
xj = axj-1 - byj-1, yj = bxj-1 + ayj.1, j = 2, ., N, (12)
где x1 = cosAa; y1 = sinAa; Aa = (n/2)/N; a = x1 - y1 ; b = 2x1 y1.
Коэффициенты ck чебышевского разложения (10) обладают замечательным свойством - по мере увеличения степени k они монотонно убывают. Назначаемую априори в формуле (10) степень полинома N следует рассматривать как верхнюю грань, заведомо завышенную, чем это необходимо. Это позволяет подобрать необходимую степень n разложения (10) для каждой функции A0(t), e(t), o(t) в соответствии с заданной погрешностью аппроксимации в. Усечение ряда (10) до степени n < N - 1 для заданной погрешности в производится по условию
N-1
sup Z abs(ck) < в. (13)
k=n
Теперь ряд (10) можно переписать так:
f (x) = Z ckTk(x). (14)
k=0
Для вычисления скорости изменения функции f(x) = df/dx эффективно использовать, как рекомендуется в [3], формулу
f(x)=(2/At)Z kckUk-1 (x) (15)
k=1
66
Геодезия
и рекуррентное соотношение для полиномов Чебышева II рода
Uo(x) = 1; Ui(x) = 2x; Uk+1 (x) = 2xUk(x) - иы(x); k = 1, n. (16)
Можно еще упростить вычисление матрицы направляющих косинусов W, соединяя два угла О и АО в один 0 по формуле
0 = О +Д0, (17)
где угол в можно считать аналогом истинного звездного времени на меридиане Гринвича в текущий момент времени T.
Тогда матрица вращения W из выражений (2) и (3) может быть представлена как произведение трех элементарных поворотов
W = R з (-0)R2(e)Ri (а), (18)
один из которых R3(-0) = R3(-(в +Д0)) описывает мгновенное сидерическое вращение Земли, два других поворота R2(в) и R1(a) моделируют совокупно прецессионно-нутационное движение оси вращения Земли как относительно ИСК, так и полярное движение относительно ОЗСК.
Два медленно и плавно меняющихся вращения R2(в) и R 1(а) выделим из формулы (18) в отдельную матрицу Q.
Q = R 2 (fi)Ri (а). (19)
Окончательно матрица W преобразования ОЗСК в ИСК запишется теперь
так:
W = R3 (-0)Q. (20)
Поскольку кардановые углы fi(t), a(t) малы на интервале [Tn, Tk], то из матрицы Q можно выделить единичную матрицу E и матрицу ДQ (описывающую совокупно колебания оси и неравномерность вращения Земли) следующим образом:
Q = E +ДQ, (21)
~-p2l 2 pG -p "
ДQ = 0 -G2l 2 G . (22)
_ в -G - (P2 +G2)/2
Тогда формулы прямого и обратного преобразования произвольного вектора из общеземной системы координат x в инерциальную у на основании (2), (20)-(22) принимают вид:
у = R з (-0)(E + ДQ)x; x = (E - AQ)R3 (0)y. (23)
67
Геодезия
Прямое и обратное преобразование координат вектора скорости из ОЗСК в ИСК производится на основании (23) по формулам:
У' = Rs (-6)(E +AQ)x' - Я3'(-в)(Е +AQ)x + R3 (-O)AQ'x; (24)
x' = (E - AQ)R3(0)y' + (E - AQ)R3 '(O)y - AQ'R3(O)y, (25)
- sin# - cos# 0 ~-№ P'a + Pa' -e "
R3'(-O) = O' cos# - sin# 0 , AQ' = 0 - aa' a'
_ в -a' -pp'-aa'
0 i о о
Выводы
1. Исходная классическая модель вращения Земли [4], представляемая формулами (1) и (2) на больших интервалах времени (более сотни лет), преобразована для ограниченных интервалов времени [Tn, Tk ] к более простому виду - формулы (20)-(25).
2. Тригонометрические члены (более двух тысяч), содержащиеся в исходной модели формулы (1), в новой модели формул (19), (20) почти полностью исключены (остались две функции sinQ и cosO).
3. Новая модель, предложенная в формулах (19), (20), адаптируясь к заданной погрешности в аппроксимации вращения Земли (и, следовательно, к точности преобразования координат) путем усечения равномерно сходящегося ряда Чебышева, является алгоритмически настраиваемой моделью по условию выражения (13).
4. С помощью новой модели проще преобразуются координаты векторов положения по формуле (23) и скорости по формулам (24), (25), как космического аппарата, так и наблюдателя, связанного с вращающейся Землей, из инерциальной системы отсчета в общеземную систему координат, и наоборот.
5. Новая модель вращения Земли, увеличивая быстродействие алгоритма численного интегрирования дифференциальных уравнений движения космических аппаратов практически без снижения точности, повышает эффективность динамического метода решения задач космической геодезии [6].
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Брумберг, В.А. Релятивистские системы координат и шкалы времени [Текст] /
В.А. Брумберг // Труды Института прикладной астрономии РАН. - Санкт-Петербург. -2004. - вып. 10. - С. 44-61.
2. Лукашова, М.В. Небесное эфемеридное начало (CEO) [Текст] / М.В. Лукашова, М.Л. Свешников // Труды Института прикладной астрономии РАН. - Санкт-Петербург. -2004. - вып. 10. - С. 186-206.
3. Глебова, Н.И. Интерполирование табличных данных [Текст] / Н.И. Глебова // Труды Института прикладной астрономии РАН. - Санкт-Петербург. - 2004. - вып. 10. - С. 39-43.
68
Геодезия
4. IERS Conventions 2003 (ed. McCarthy D. D. and G. Petit) [Текст] / IERS Technical Note № 32. - Frankfurt am Main. - 2004. - 127 c.
5. Каула, У. Спутниковая геодезия. Теоретические основы [Текст] / У. Каула. - М.: Мир, 1976. - 172 с.
6. Сурнин, Ю.В. Программный комплекс «Орбита-СГГА-2» для решения задач космической геодезии динамическим методом [Текст] / Ю.В. Сурнин, В.А. Ащеулов, С.В. Куже-лев, Е.В. Михайлович, Н.К. Шендрик // Геодезия и картография. - 2008. - № 2. - С. 14-19.
7. Марчук, Г.И. Методы вычислительной математики [Текст] / Г.И. Марчук. - Новосибирск: Наука, 1973. - 352 с.
8. Демидович, Б.П. Численные методы анализа [Текст] / Б.П. Демидович, И.А. Марон, Э.Э. Шувалова. - М.: Физ.-мат. лит., 1963. - 400 с.
Получено 25.06.2010
© Ю.В. Сурнин, 2010
69
