Примеры алгоритмов расчета координат и скорости нка по данным эфемерид и альманаха

Ниже даны примеры алгоритмов расчета координат и составляющих скорости НКА ГЛОНАСС на текущий момент времени по данным эфемерид и альманаха системы.

П.3.1 Пример алгоритма пересчета эфемерид НКА на текущий момент времени

Пересчет эфемерид с момента времени t_э на моменты t_i измерения навигационных параметров (½t _i ½= ½t _i - t_э ½ < 15 мин ) проводится методом численного интегрирования дифференциальных уравнений движений КА, в правых частях которых учитываются ускорения, определяемые константой гравитационного поля Земли m, второй зональной гармоникой С₂₀ , характеризующей полярное сжатие Земли, а также ускорения лунно-солнечных гравитационных возмущений.

Уравнения движения интегрируются в прямоугольной абсолютной геоцентрической системе координат OX_aY_aZ_a, связанной с текущими экватором и точкой весеннего равноденствия, методом Рунге-Кутта четвертого порядка и имеют вид:

dxa/dt = Vxa dya/dt = Vya dza/dt = Vza _ _ _ _ _ dVxa/dt = - m * Xa + 3/2 * C20 * m * Xa *r2 * (1 - 5 * Za2 ) + Jxaл + Jxaс, _ _ _ _ _ dVya/dt = - m * Ya + 3/2 * C20 * m * Ya *r2 * (1 - 5 * Za2 ) + Jyaл + Jyaс, _ _ _ _ _ dVza/dt = - m * Za + 3/2 * C20 * m * Za *r2 * (1 - 5 * Za2 ) + Jzaл + Jzaс,

    ( 1 )   

_ _ _ _

Здесь m = m / r² , X_a = x_a/r , Y_a = y_a / r , Z_a = z_a/r , r = a_e / r ,

_____________

r = Ö X_a² + Y_a² + Z_a² .

Jxaс,Jyaс,Jzaс	-	ускорения от солнечных гравитационных возмущений ;
Jxaл,Jyaл,Jzaл	-	ускорения от лунных гравитационных возмущений;
ae	-	экваториальный радиус Земли, равный 6378,136 км [ПЗ-90, справочный документ, КНИЦ, 1998];
m	-	константа гравитационного поля Земли ( 398600,44 км3/c2 ) [ПЗ-90];
С20	-	коэффициент при второй зональной гармонике разложения гравита-ционного поля Земли в ряд по сферическим функциям, равный минус 1082,63*10-6 ( 20 = * 20, где 20 - нормализованное значение гармонического коэффициента при второй зональной гармонике, равное минус 484,165*10-6 [ПЗ-90]).

Ускорения от лунных и солнечных гравитационных возмущений вычисляются по формулам:

_ _ Jxak = mk [ ( xк - Xak ) / D3к - xк ] , _ _ Jyak = mk [ ( hк - Yak ) / D3к - hк ] , _ _ Jzak = mk [ ( zk - Zak ) / D3к - zк ] ,

   ( 2 )  

_ _ _ _

где: m _k = m _k / r²_k , X_ak = X_a / r_k , Y_ak = Y_a / r_k , Z_ak = Z_a / r_k ,

_ _ _

D²_k = (x_k - X_ak )² + (h_k - Y_ak )² + (z _k - Z_ak )²,

k	-	индекс возмущающего тела, k = л для Луны и к = с для Солнца;
xk , hk , zk , rk	-	направляющие косинусы и радиус-вектор возмущающих тел в системе OXaYaZa на момент tэ
mл	-	константа гравитационного поля Луны, равная 4902,835 км3/с2
mс	-	константа гравитационного поля Солнца, равная 0,1325263 * 1012 км/с2

Входящие в (2) величины x_k , h_k , z_k , r_k вычисляются один раз (на момент времени t_э ) на весь интервал размножения ( ± 15 мин) по формулам [Дубошин Г.Н., Небесная механика: Основные задачи и методы; М.: Наука, 1975; Абалакин В.К., Основы эфемеридной астрономии, М.: Наука, 1979]:

xл = sin(л + Г’) x11 + cos(л + Г‘) x12 , hл = sin(л + Г’) h11 + cos(л + Г‘) h12 , zл = sin(л + Г’) z11 + cos(л + Г‘) z12 , xсэ = cos c * cos wc - sin c * sin wc , hc = (sin c * cos wc + cos c * sin wc ) cos e , zc = (sin c * cos wc + cos c * sin wc ) sin e , rk = ak ( 1 - ek * cos Ek ), ( k = л, с )

    ( 3 )   

где E_k = g_k + e_k * sin E_k ,

______

sin _k =Ö`1- e_k² * sin E_k * ( 1 - e_k * cos E_k )^-1 ,

cos _k = ( cos E_k - e_k ) * ( 1 - e_k * cos E_k )^-1 ,

x₁₁ = sin W_л * cos W_л * ( 1 - cos i_л ) ,

x₁₂ = 1 - sin² W_л * ( 1 - cos i_л ) ,

h₁₁ = x^* * cos e - z^* * sin e ,

h₁₂ = x₁₁ * cos e + h^* * sin e ,

z₁₁ = x^* * sin e + z^* * cos e ,

z₁₂ = x₁₁ * sine - h^* * cose ,

x^* = 1 - cos² W_л ( 1 - cos i_л )

h^* = sin W_л * sin i_л ,

z^* = cos W_л * sin i_л ,

g_k = g_ok + g_1k * T,

W_л = W_ол + W_1л * T,

Г^‘ = Г’₀ + Г’₁ * Т,

Т = ( 27392,375 + Sдн + t_э / 86400 ) / 36525,

где:

ал	-	большая полуось орбиты Луны, равная 3,84385243*105 км;
ас	-	большая полуось “орбиты” Солнца, равная 1,49598*108 км
ел	-	эксцентриситет лунной орбиты, равный 0,054900489;
ес	-	эксцентриситет солнечной орбиты, равный 0,016719
iл	-	наклонение орбиты Луны к плоскости эклиптики, равное 508'43'',4;
e	-	средний наклон эклиптики к экватору, равный 2326'33'';

g_ол = -63 53' 43'',41

g_1л = 477198 50' 56'',79

W_ол = 259 10' 59'',79

W_1л = -1934 08' 31'',23

Г’_о = -334 19' 46'',40

Г’₁ = 4069 02' 02'',52

w_c = 281 13' 15'',00 + 6189'',03 * T;

g_ол = 358 28' 33'',04;

g_ол = 129596579'',10;

Т - время от основной эпохи 1900, янв. 05 (GMT) до момента задания эфемерид t_э в юлианских столетиях по 36525 эфемеридных суток;

27392,375 - число дней от основной эпохи 1900, янв 05 (GMT) до эпохи 1975, янв. 0. (МДВ) с учетом трех часов при пересчете московского декретного времени (МДВ) t_э в гринвичское (GMT);

Sдн - сумма дней от 0 ч эпохи 1975, янв.0 (МДВ) до 0ч текущей даты (МДВ), к которой относится время t_э ( отсчет начала дат по московскому времени).

Начальными условиями для интегрирования системы (1) являются гринвичские координаты X(t_э), Y(t_э), Z(t_э) и составляющие вектора скорости Vx(t_э) , Vy(t_э), Vz(t_э), содержащиеся в навигационном кадре, которые пересчитываются из гринвичской системы координат OXYZ (ПЗ-90) в абсолютную OX_aY_aZ_a по формулам :

X_a(t_э) = X(t_э) * cosS - Y(t_э) * sinS ,

Y_a(t_э) = Y(t_э) * sinS + Y(t_э) * cosS,

Z_a(t_э) = Z(t_э),

Vx_a(t_э) = Vx(t_э) * cosS - Vy(t_э) * sinS - w_з* Y_a(t_э),

Vy_a(t_э) = Vx(t_э) * sinS + Vy(t_э) * cosS + w_з* X_a(t_э),

Vz_a(t_э) = V_z(t_э) ,

S = s + w_з ( t - 3^h ).

Здесь

w_з- угловая скорость вращения Земли, равная 0,7292115 * 10^-4 с^-1;

s - истинное звездное время в гринвичскую полночь даты задания эфемерид t_э.

Примечания.

1. Ускорения Jx_aс, Jx_aл, Jy_aс, Jy_aл, Jz_aс, Jz_aл в (1) могут быть приняты постоянными и вычисляться один раз на момент t_э по формулам (2) или исключены из (1) с последующим добавлением к результатам интегрирования поправок

DX = ( Jх_aл + Jх_aс ) * t²/2 , DY = ( Jy_aл + Jy_aс ) * t²/2 , DZ =( Jz_aл ₊ Jz_aс ) t²/2 ,

DVx = ( Jх_aл + Jх_aс ) * t , DVy = ( Jy_aл + Jy_aс ) * t , DVz =( Jz_aл ₊ Jz_aс ) t ,

где  = t_i - t_э.

2. Направляющие косинусы x_к , h_к , z_k могут вычисляться по формулам (3) или передаваться извне.

3. Начало гринвичской ( правой ) системы координат - в центре масс Земли, ось OZ направлена по оси вращения Земли к северному полюсу, а ось OX - в точку пересечения гринвичского меридиана с плоскостью экватора.

4. Если при интегрировании системы (1) исключить лунно-солнечные ускорения (2), а их учет производить добавлением к результатам интегрирования поправок

DX = ( Jх_aл + Jх_aс ) * t²/2 , DY = ( Jy_aл + Jy_aс ) * t²/2 , DZ =( Jz_aл ₊ Jz_aс ) t²/2 ,

DVx = ( Jх_aл + Jх_aс ) * t , DVy = ( Jy_aл + Jy_aс ) * t , DVz =( Jz_aл ₊ Jz_aс ) t ,

то возникающее при этом увеличение ошибок размножения эфемерид не превышает 10 % .

Здесь (Jх_aл + Jх_aс), (Jy_aл + Jy_aс) , (Jz_aл ₊ Jz_aс) - проекции лунно-солнечных ускорений на оси системы OX_aY_aZ_a на момент задания эфемерид t_э вычисляются по формулам (2).

5. Для расчета эфемерид КА на моменты навигационных измерений t_j можно использовать проекции лунно-солнечных ускорений X(t_э), Y(t_э), Z(t_э) на оси гринвичской геоцентрической системы координат, которые передаются в составе навигационного кадра. Перед интегрированием системы дифференциальных уравнений (1) эти ускорения должны быть переведены в прямоугольную абсолютную геоцентрическую систему координат OX_aY_aZ_a по формулам:

(Jх_aл + Jх_aс) = X(t_э) * cos S - Y(t_э) * sin S ,

(Jy_aл + Jy_aс) = X(t_э) * sin S + Y(t_э) * cos S ,

(Jz_aл ₊ Jz_aс) = Z(t_э)

В таблице приведены величины точности размножения эфемерид (в метрах)

Шаг интегрирования (мин.)	Интервал интегрирования
	5 мин	10 мин	15 мин
1	0.42	0.56	0.77
2.5	0.42	0.56	0.77
5	0.45	0.61	0.83
7.5	-	-	1.21

П.3.2 Алгоритм расчета параметров движения спутников по данным альманаха

Алгоритм расчета параметров движения НКА ГЛОНАСС по данным альманаха системы (АС) используется потребителем при выборе оптимального созвездия, расчете целеуказаний для вхождения в связь с выбранным НКА. Назначение алгоритма - расчет координат и составляющих вектора скорости НКА на моменты t_i вхождения потребителем в связь с НКА.

.

П.3.2.1 Состав данных, образующих АС

АС содержит набор параметров орбит НКА системы ГЛОНАСС, заданных для каждого НКА на момент прохождения им первого (внутри суток с номером N_Aj) восходящего узла орбиты t_lj.

Набор параметров орбит каждого НКА содержит:

NAj	-	календарный номер суток внутри четырехлетнего периода от начала ближайшего високосного года, к которым относятся данные АС для j-го НКА ;
lj	-	гринвичская долгота восходящего узла орбиты j-го НКА момент tlj (радианы);
tlj	-	московское декретное время прохождения j-м НКА восходящего узла орбиты, ближайшее к началу суток с номером NAj (секунды);
Dij	-	поправка к среднему значению наклонения орбиты j-го НКА на момент tlj ;
DTj	-	поправка к среднему значению драконического периода обращения j-го NKA ;
DТj	-	скорость изменения периода обращения j-го НКА ;
ej	-	эксцентриситет орбиты j-го НКА на момент времени tlj ;
wj	-	аргумент перигея орбиты j-го НКА на момент времени tlj (радианы).

Здесь l - индекс принадлежности параметров АС ко времени прохождения восходящего узла орбиты t_lj , а j - номер НКА (j=1,......,24). В дальнейшем индекс j опущен.

Средние значения наклонения плоскости орбиты НКА системы ГЛОНАСС i_ср и драконического периода обращения Т_ср составляют 63⁰ и 43200 с, соответственно.

Набор параметров орбиты для каждого НКА задан в гринвичской геоцентрической системе координат OXYZ, “замороженной” в момент t.

Начало системы совмещено с центром масс Земли. Ось Z направлена к северному полюсу на среднюю эпоху 1900-1905 г.г., плоскость XOZ при этом параллельна среднему гринвичскому меридиану и определяет положение нуль-пункта системы счета долгот, ось OY дополняет систему до правой.

П.3.2.2 Алгоритм расчета

Расчет координат и составляющих вектора скорости НКА по данным альманаха системы ГЛОНАСС в абсолютной геоцентрической системе координат OX_aY_aZ_a (начало системы координат и направление оси OZ_a совпадает с началом системы координат OXYZ и направлением оси OZ, плоскость XOZ отстоит от плоскости X_aOZ_a на величину истинного звездного времени S, а ось OY_a дополняет систему до правой) на заданный момент времени t_i (московское декретное время суток с номером N_o внутри четырехлетнего периода) проводится в два этапа.

Сначала с помощью величин DТ, DТ и l рассчитываются момент прохождения восходящего узла орбиты t_k на витке с номером K, к которому принадлежит заданный момент времени t_i (t_i - t_k < T_ср + DТ ), и долгота восходящего узла l_k на этом витке. Остальные параметры принимаются постоянными и равными тем, которые содержатся в навигационном кадре.

Затем оскулирующие элементы пересчитываются с момента t_k по аналитическим формулам на момент t_i . При этом учитываются вековые и периодические возмущения в элементах орбиты НКА от второй зональной гармоники C₂₀ в разложении геопотенциала, характеризующей полярное сжатие Земли.

Полученные на момент t_i оскулирующие элементы переводятся в кинематические параметры. Последовательность проведения расчета и используемые рабочие формулы приведены ниже.

1) Методом последовательного приближения находится большая полуось орбиты а:

а⁽ⁿ⁺¹⁾= [m * (T / 2p )² ]^1/3,

Т_оск⁽ⁿ⁺¹⁾ = Т_др *{ 1+3/2 C₂₀ (а_е/p⁽ⁿ⁾)² *

*[ (2 - 5/2sin² i) * (1-e²)^3/2 / (1 + e*cos w )² + (1+ e* cos )³ / (1-e²)]}^-1,

p⁽ⁿ⁾ = a⁽ⁿ⁾ * ( 1- e² ), n = 0, 1, 2...,

где  = -w, i = i_ср + Di и Т_др = Т_ср + DТ.

За начальное приближение принимается а⁽⁰⁾= [m * (Т_др / 2p )² ]^1/3,

Приближение заканчивается при выполнении условия ½а^(n-1) - a⁽ⁿ⁾½ < 10^-3 км.

Для этого обычно достаточно двух-трех итераций.

2) Рассчитываются момент прохождения восходящего узла орбиты t_k на витке, к которому принадлежит момент t_i, и долгота восходящего узла на этом витке l_k:

_

t_lk = [ t_lk ] _{mod 86400}

_

t_lk = t_l + T_др * W + DТ * W²

W_k = t^* / T_др, W - целая часть W_k,

t^* = t_i - t_l + 86400 * (N_o - N_A) ,

l_k = l + (W - w₃) * ( W * T_др + DТ * W² ),

W = 3/2 * C₂₀ * n * ( a_e / a )² * cos i * (1-e²)^-2,

n = 2 p / Т_др , W = l_k + S , S = S₀ + w₃ (t_lk - 10800 ).

Здесь:

С20	-	коэффициент при второй зональной гармонике разложения геопотенциала в ряд по сферическим функциям, равный -1082,63 * 10-6,
ае	-	экваториальный радиус Земли, равный 6378,136 км,
S0	-	истинное звездное время на гринвичскую полночь даты N0 , к которой относится время ti ,
w3	-	угловая скорость вращения Земли, равная 0,7392115 * 10-4 с-1,
	-	константа гравитационного поля Земли, равная 398600,44 км3 / с2.

3) Вычисляются константы интегрирования на момент t_lk :



dа(m) / а	=	2*J* (ae / a)2*(1-3/ 2sin2 i)*( l * cos L ) + J*( ae / a )2*sin i * * (1/2 * h * sin L - 1/2* l * cos L + cos 2l + 7/2 * l * cos 3L + 7/2*h*sin3L)
dh(m)	=	J*(ae / a)2*(1-3/ 2sin2 i)*[ l*n* t + sin 3L + 3/2 * l * sin 2L - 3/2*h*cos2L] - - 1/4*J*(ae / a )2*sin2 i*[ sin L - 7/3*sin3L + 5*l*sin 2L -17/2 * l * sin 4L + + 17/2*h*cos4L + h*cos2L] + J*(ae / a )2*cos2 i*(l*n* t- 1/2 * l * sin 2L)
dl(m)	=	J*(ae / a)2*(1-3/2 sin2 i)*[-h*n*t + cos L + 3/2 * l * cos 2L + 3/2*h*sin 2L] - - 1/4*J*(ae / a )2*sin2 i*[-cos L - 7/3*cos3L - 5*h*sin 2L - 17/2 * l * cos 4L - - 17/2*h*sin4L + l*cos2L] + J*( ae / a )2*cos2 i*(-h*n*t + 1/2 *h * sin 2L)	(1)
dW(m)	=	J*(ae / a)2*cos i*( n*t + 7/2 * l * sin L - 7/2*h*cosL -1/2sin2L - 7/6*sin3L + + 7/6*h*cos 3L)
di(m)	=	1/2*J*(ae / a)2*sin i * cos i * ( -1 *cos L +h*sin L +cos2L +7/3*l*cos 3L + + 7/3*h*sin 3L)
dL(m)	=	2*J*(ae / a)2*(1 - 3/2*sin i )*( n*t + 7/4 * l * sin L - 7/4 * h * cos L) + + 3*J*(ae / a)2*sin i * (-7/24 * h * cos L - 7/24 * l * sin L -49/72* h *cos 3L + + 49/72* l *sin 3L + 1/4 * sin 2L) + J * (ae / a)2* cos i * ( n*t + 7/2 * l * sin L - - 5/2 h cos L - 1/2 * sin 2L - 7/6 * l*sin 3L + 7/6 * h*cos 3L)



где L = M + , M = E -  * sin E, tg(E/2) = [(1 - ) / (1 + )]^1/2 *tg(/2),

h =  * sin , l =  * cos , m = 1,

 = 0, J = 3/2 * C₂₀, a = a⁽ⁿ⁾ (из пункта 1).

4) Вычисляются поправки к элементам орбиты НКА на момент времени t_i за счет влияния второй зональной гармоники С₂₀:

	dа = dа(2) - dа(1)
	dh = dh(2) - dh(1)
	dl = dl(2) - dl(1)
	dW = dW(2) - dW(1)
	di = di(2) - di(1)
	dL* = dL(2) - dL(1)

Величины dа⁽²⁾, dh⁽²⁾, dl⁽²⁾, dW⁽²⁾, di⁽²⁾, dL⁽²⁾ вычисляются для  = t_i - t_lk и m =2 по формулам (1) при L = M +  + n * .

5. Вычисляются возмущенные элементы орбиты НКА на момент времени t_i:

	ai = a + dа,
	hi = h + dh,
	li = l + dl,
	i = (hi * hi + li * li)1/2,
	i=	{	arctg (hi / li ), если i  0 и li  0 0, если i = 0,  / 2, если i  0 и li = i, - / 2, если i  0 и li = -i,
	Wi = W + dW
	ii = i + di
	Mi = L* - i , L* = M +  + n * (ti - tlk) + dL*.

Здесь i - индекс принадлежности ко времени t_i,

6. Вычисляются координаты и составляющие вектора скорости НКА в системе координат OX_aY_aZ_a на момент времени t_i :

	Ei(n) = Mi + i * sin Ei(n-1), Ei(0) = Mi,  Ei(n) - Ei(n-1)   10-8,
	tg(i/2) = [(1 + i) / (1 - i)]1/2 *tg(Ei(n)/2), ui = i + i,
	ri = ai * ( 1 - i * cos Ei(n)),
	Vri = (  / ai)1/2 * (i - sin i) / (1 - i * i)-1,
	Vui = (  / ai)1/2 * (1 + i * cos i) / (1 - i * i)-1,

Xi	=	ri * ( cos ui * cos i - sin ui * sin i * cos ii),
Yi	=	ri * ( cos ui * sin i + sin ui * cos i * cos ii),
Zi	=	ri * sin ui * sin ii,
Vxi	=	Vri * ( cos ui * cos i - sin ui * sin i * cos ii) - -Vui * ( sin ui * cos i + cos ui * sin i * cos ii),
Vyi	=	Vri * ( cos ui * sin i + sin ui * cos i * cos ii) - -Vui * ( sin ui * sin i - cos ui * cos i * cos ii),
Vzi	=	Vri * sin ui * sin ii + Vui * cos ui * sin ii.

По всем вопросам о состоянии и использовании системы ГЛОНАСС Вы можете обращаться в Координационный научно-информационный центр.

Россия, 117279, г. Москва, а/я 14

Телефон (факс): (095) 333-81-33

e-mail: sfcsic@mx.iki.rssi.ru

Internet: http://www.rssi.ru/SFCSIC/SFCSIC_main.html

© 1998 КООРДИНАЦИОННЫЙ НАУЧНО-ИНФОРМАЦИОННЫЙ ЦЕНТР

39