CC BY
127
УДК 528.33 Ю.В. Сурнин СГГА, Новосибирск
АППРОКСИМАЦИЯ МАТЕМАТИЧЕСКОЙ МОДЕЛИ ВРАЩЕНИЯ ЗЕМЛИ УГЛАМИ КАРДАНО И ПОЛИНОМАМИ ЧЕБЫШЕВА
Рассматривается замена двенадцати углов, описывающих в классической модели эффекты прецессии, нутации, полярного движения и неравномерности вращения Земли, тремя кардановыми углами, которые аппроксимируются полиномами Чебышева. В результате более двух тысяч тригонометрических и полиномиальных членов заменяются небольшим числом чебышевских многочленов. Построенная таким образом модель вращения Земли повышает быстродействие алгоритма, практически не снижает точность и потому более эффективна для решения задач космической геодезии динамическим методом.
Yu.V. Surnin SSGA, Novosibirsk
CARDAN ANGLE-AND CHEBYSHEV POLYNOMIAL APPROXIMATION OF THE EARTH ROTATION SIMULATOR
The author considers substitution of the nine angles describing (in classical model) the effects of precession, nutations, polar motion and non-uniformity of the Earth rotation for three cardan angles approximated by Chebyshev polynomials. As a result more than a thousand of trigonometric and polynomial members are replaced by a small number of Chebyshev multinomials. The Earth rotation model thus built is more efficient as the algorithmic speed increases without considerable degradation of the model accuracy.
Математическая модель вращения Земли обычно представляется в виде произведения элементарных матриц поворота [1, 2, 4], учитывающих движение оси мира относительно инерциального пространства, собственно вращение Земли, неравномерность ее вращения и полярное движение мгновенной оси в теле Земли. Чтобы уменьшить влияние погрешностей модели вращения Земли на орбиты космических аппаратов, целесообразно в качестве инерциальной системы координат (ИСК) - точнее квазиинерциальной - выбрать небесную систему отсчета [4], определяемую истинным равноденствием и истинным экватором фиксированной эпохи Tm, лежащей, примерно, в середине некоторого интервала [Tn, Tk]. Например, Tn, Tk - эпохи начала и конца какой-либо наблюдательной компании. Тогда матрица W направляющих косинусов осей общеземной системы координат (ОЗСК) относительно осей ИСК в рамках ньютоновой механики может быть описана двенадцатью элементарными поворотами
W= Ri (-Som-ASm)R3(-^^m)Rl(S0m)R3(ZAm)R2(-OAm)R3(ZAm)Rl (-ео)Яз(А ^JRjfSo+As)^ ^R3(-S)R1(xp)R2(yp). (1)
В формуле Аут, Ау, и Ает, Ае - компоненты астрономической нутации в долготе и наклоне соответственно в эпохи Тт и Т; Слт, вАт, 2Ат -экваториальные параметры прецессии на интервале времени [Тт, Т]; е0, и е0т - средние наклоны эклиптики к экватору в эпохи Т и Тт; 5 - истинное звездное время на меридиане Гринвича на момент времени Т; хр, ур -параметры полярного движения оси вращения в теле Земли в эпоху Т. Параметр неравномерности вращения Земли АиТ1(Т) включен в угол 5.
Если вместо ньтоновой механики для модели вращения Земли использовать общую теорию относительности [1, 4], то дальнейшие преобразования матрицы W, о которых пойдет речь в статье, принципиально не изменяются. Меняется только сама исходная матрица W в выражении
(2)
которое осуществляет связь прямоугольных координат произвольного вектора х={х1, х2, х3},заданного в ОЗСК, с его координатами у={у1, у2, у3} в ИСК.
Современная аналитическая модель вращения Земли [1, 2, 4] содержит более тысячи полиномиальных и тригонометрических членов. Синусы и косинусы значительно «тормозят» процесс численного интегрирования уравнений движения искусственных спутников Земли в динамическом методе космической геодезии [6]. С целью увеличения скорости расчетов предлагается преобразовать традиционную формулу (1) к виду
W=Rз(Ш■ (3)
В правой части (3) выделен главный член R3(-в), содержащий некоторый аналог гринвичского звездного угла в [1, 4]. Угол в вводится, как функция земного времени ТТ (в [2] используется всемирное время ЦТ1), практически реализуемого атомными стандартами частоты, по формулам:
в= вт + ю(Т-Тт), От=во+ы(Тт-То) (4)
во=2п0,7790572732640радиан, т=2ж(1+^)радиан/сутки,
¡л=0,00273 7811911354480,
где т - номинальная угловая скорость вращения Земли, То=2451545.0 -стандартная эпоха ^000.0.
Остальная часть в (3) - матрица Ц - описывает небольшие вековое движение и колебания оси вращения Земли на интервале времени [Ти, Т] выражением
Q=Rз(в)W. (5)
С одной стороны, матрицу направляющих косинусов Ц можно представить как матрицу
дг=\1’/ к’], (6)
строками которой являются орты I, у , к осей ИСК в проекциях на орты ОЗСК
1={1, 0, 0},у={0, 1, 0}, к={0, 0, 1}.
С другой стороны, эту же матрицу Ц можно представить в виде произведения трех матриц элементарных поворотов либо на углы Эйлера, либо на углы Кардано. При малых колебаниях экватора (по наклону) возникает квазилинейная зависимость двух углов Эйлера (прецессии и чистого вращения), примыкающих друг к другу по линии узлов. Отдавая предпочтение углам Кардано, правую часть (6) представим как
Q=Rз(-Ав)R2(P)Rl(a), (7)
где кардановые углы можно интерпретировать следующим образом: Ав -аналог нутации по прямому восхождению, в - аналог нутации по склонению и а- аналог нутации по наклону. Углы Кардано Ав, в, & можно вычислить через исходные орты I, у , к , получаемые как столбцы матрицы Ц, по формулам (1), (5) и (6). Если ввести три вспомогательных орта I, ] , к равенствами:
]=кМ/\ к’\, к=1х], 1=]^к, (8)
то в области определения [-п/2, п/2] углы Кардано найдутся по формулам:
tgАв=-ij /И , tgв=k’^i/k ^ , tgo=kj /№ . (9)
Следует обратить внимание на то, что эмпирические углы хр, ур и (1 +^)А иТ1 даются не аналитическими выражениями, как это имеет место для всех остальных углов в равенстве (1), а публикуются в виде таблицы по аргументу Т. Поэтому табличное задание трех функций хр(Т), ур(Т), АиТ1(Т) необходимо преобразовать в аналитический вид, используя, например, аппроксимацию кубическими сплайнами [7].
Теперь малые углы Кардано Ав(?), в^), в@) можно рассматривать как аналитические функции времени t, даваемые формулами (1), (5), (6), (8), (9) и аппроксимировать их на интервале [Ти, Тк] полиномами Чебышева, используя чебышевский альтернанс и равномерную сходимость ряда [8].
Пусть /(х) любой из углов Кардано Ав(0, в^), &(^, t=T-Tm. Тогда разложение функции /(х) по полиномам Чебышева до степени N-1 можно записать так
N-1
/0) = ЕсЛ(х)’ (10)
к=0
где х=(Ы0)Ш, ^=^+^)/2, Аt=(tk-tn)2, t=T-Tm , -1< х < 1, Т0(х)=1, Т1(х)=х, Ти+1(х)=2хТи(х)-Ты(х), k=1,...,N-1.
Коэффициенты ck разложения (10) находятся по формулам:
j~N
ск=[(2-д0кУЩ £ f(Xj)Tk(Xj), к=0, 1, 2,..., N-1, док 1, если к 0 иначе 0, (11)
где функции f(xj) вычисляются с помощью равенств (1), (5), (6), (8), (9) в корнях Xj полинома Чебышева TN степени N. Заметим, что cN=0, поэтому в ряде (10) предел суммирования равен N-1. Для всех трех функций A0(t), fi(t), a(t) узлы Xj вычисляются один раз (с минимальным обращением к функциям cos и sin) по рекуррентным формулам:
j=axj-i-byj-i, yJ=bxJ-i+ayJ-i, j=2,...,N, (12)
2 2
где x1=cosAa, y1=sinAa, Aa=(n/2)/N, a=x1 -y1 , b=2x1y1.
Для вычисления скорости изменения функции f(x)=df/dx эффективно использовать, как рекомендуется в [2], формулу
АЧ
f(x)=(2/At)Jj кскиы(х), (13)
Ы
и рекуррентное соотношение для полиномов Чебышева II рода
По(х)=1, и1(х)=2х, ик+1(х)=2хЩх) - П^х), k=1,..., N-1. (14)
Коэффициенты ск чебышевского разложения (10) обладают
замечательным свойством - по мере увеличения степени к они монотонно убывают. Назначаемую априори степень полинома N следует рассматривать как максимальную, т.е.: заведомо завышенную, чем это необходимо. Это позволяет подобрать необходимую степень п разложения (10) в соответствии заданной погрешностью аппроксимации в. Усечение ряда (10) до степени п<Ы-1 для заданной погрешности в производится по условию
N-1
supY, abs(cjj<e. (15)
Поскольку кардановые углы А0^), ¡3^), малы на интервале [Тп, Тк], то
из матрицы Ц можно выделить единичную матрицу Е и матрицу АЦ (описывающую совокупно колебания оси и неравномерность вращения Земли) следующим образом:
ЦЕ+АЦ, (16)
(А<92 + р2)/2 А0 в
- А0 + р<7
-(А#2 +о2)/2 - о
-Р-сгА0 О-РА0 -(в V о 2)/2
k=
n
Тогда формулы прямого и обратного преобразования произвольного вектора из общеземной системы координат x в инерциальную у принимают вид:
y=R3(-0)(E+AQ)x, x=(E-AQ)R3(6)y. (19)
С помощью выражений (19) упрощается прямое и обратное преобразование координат вектора скорости из ОЗСК в ИСК или наоборот. Дифференцируя по времени t выражения (19), получаем
У’= Rs(-e)(E+AQ)xf - 0'R3' (-d)(E +AQ)x + R3(-d)AQ'x, (20)
x'=(E-AQ)R3(0)y' + в' (E-AQ)Rs' (в)у-AQ'R3(0)y, (21)
где x'=dx/dt, y = dy/dt, ß'=dß/dt, AQ'=dAQ/dt, R3' (0)=dR3(0)/dt.
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Брумберг В.А. Релятивистские системы координат и шкалы времени [Текст] / В.А. Брумберг // Труды Института прикладной астрономии РАН. - Санкт-Петербург. -2004.- вып. 10. - С. 44-61.
2. Лукашова М.В., Свешников М.Л. Небесное эфемеридное начало (CEO) [Текст] / М.В. Лукашова, М.Л. Свешников // Там же. - С. 186-206.
3. Глебова Н.И. Интерполирование табличных данных [Текст] / Р.И. Глебова // Там же.- С. 39-43.
4. IERS Convensions 2000 (ed. McCarthy D. D.) [Текст] /IERS Technical Note № 32. -Obs. de Paris. - 2003. - 135 c.
5. Каула У. Спутниковая геодезия. Теоретические основы [Текст]/ Каула У. Изд-во «Мир». - М. - 1976. - 172 с.
6. Сурнин Ю.В. и др. Программный комплекс «Орбита - СГГА-2» для решения задач космической геодезии динамическим методом [Текст] / Ю.В. Сурнин, В.А. Ащеулов, С.В. Кужелев, Е. В. Михайлович, Н. К. Шендрик. // Геодезия и картография. -2008. - № 2. - С. 14-19.
7. Марчук Г.И. Методы вычислительной математики [Текст] / Г.И. Марчук. Изд-во «Наука». - Сибирское отд. - Новосибирск. - 1973. - 352с.
8. Демидович Б.П. и др. Численные методы анализа [Текст] / Б.П. Демидович, И.А. Марон, Э.Э. Шувалова // Изд-во «Физ. - мат. лит.» - М. - 1963. - 400 с.
© Ю.В. Сурнин, 2010
