2013-09-12 4 views
2

Я пытаюсь рассчитать скорость полета спутников с помощью Python и pyephem. К сожалению, результат pyephems кажется неправильным.Неправильный диапазон скорости с Pyephem

После сравнения значения с расчетами, выполненными другими программами спутникового слежения, такими как GPredict или Ham Radio Deluxe, разница увеличивается до 2 км/с. Вычисленные значения для лодыжки Azemuth и Elevation - это почти та же мысль. TLE - новые, а системные часы - то же самое.

Вы видите ошибку в моем коде или у вас есть идея, что еще может вызвать ошибку?

спасибо!

Вот мой код:

import ephem 
import time 
#TLE Kepler elements 
line1 = "ESTCUBE 1"  
line2 = "1 39161U 13021C 13255.21187718 .00000558 00000-0 10331-3 0 3586" 
line3 = "2 39161 98.1264 332.9982 0009258 190.0328 170.0700 14.69100578 18774" 
satellite = ephem.readtle(line1, line2, line3) # create ephem object from tle information 
while True: 
    city = ephem.Observer() # recreate Oberserver with current time 
    city.lon, city.lat, city.elevation = '52.5186' , '13.4080' , 100 

    satellite.compute(city) 
    RangeRate = satellite.range_velocity/1000 # get RangeRate in km/sec 
    print ("RangeRate: " + str(RangeRate)) 
    time.sleep(1) 

Я записал несколько значений диапазоне скоростей от сценария и от GPRedict сделать ошибку воспроизводимо:

ESTCUBE 1    
1 39161U 13021C 13255.96108453 .00000546 00000-0 10138-3 0 3602 
2 39161 98.1264 333.7428 0009246 187.4393 172.6674 14.69101320 18883 

date: 2013-09-13 
time  pyephem-Script Gpredict   
14:07:02 -1.636   -3.204 
14:12:59 -2.154   -4.355 
14:15:15 -2.277   -4.747 
14:18:48 -2.368   -5.291 

И я добавил несколько строк, чтобы вычислить спутники высота и координаты:

elevation = satellite.elevation 
sat_latitude = satellite.sublat 
sat_longitude = satellite.sublong 

Результаты со штампом времени:

2013-09-13 14:58:13 
RangeRate: 2.15717797852 km/s 
Range: 9199834.0 
Sat Elevation: 660743.6875 
Sat_Latitude: -2:22:27.3 
Sat_Longitude: -33:15:15.4  

2013-09-13 14:58:14 
RangeRate: 2.15695092773 km/s 
Range: 9202106.0 
Sat Elevation: 660750.9375 
Sat_Latitude: -2:26:05.8 
Sat_Longitude: -33:16:01.7 

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

f_Doppler_corrected = (c0/(c0 + RangeRate))*f0 

Range Rate описывает скорость движущегося объекта на зрительной оси наблюдателя. Может быть, range_velocity - это что-то другое?

+0

какая версия Python? –

+0

Я использую Python 2.7.4 – Sebastian

+1

Чтобы помочь вам, нам нужны три вещи: точное время; положение спутника и скорость диапазона, возвращаемые этими другими инструментами, когда в то же время дается точно такой же TLE; и положение и диапазон спутника, возвращенные этим скриптом Python в тот момент. Затем мы можем, во-первых, воспроизвести ваш результат локально с использованием одного и того же сценария; затем перейдите к изучению разницы в результатах. Благодаря! –

ответ

2

Кажется, что pyephem (libastro как бэкэнд) и gpredict (предсказать) в качестве бэкэнд используют разные способы расчета скорости спутника. Я прилагаю подробный вывод сравнения для фактического контрольного наблюдения. Можно видеть, что оба выводят правильную позицию, а только gpredict выдает разумные значения range_rate. Кажется, что ошибка встречается в векторе скорости спутника. Я бы сказал, что причины из gpredict более разумны (и аналогичный код с вопросительными знаками в libastro ..), поэтому я предложу исправить в libastro, чтобы обрабатывать его, как в gpredict, однако, возможно, кто-то, кто понимает математику за ним, может добавьте к этому.

Я добавил еще один инструмент, PyPredict (также на основе прогноза), чтобы получить некоторые вычисления здесь. Однако эти значения отключены, поэтому должно быть что-то еще.

Pyephem: 3.7.5.3 
Gpredict: 1.3 
PyPredict 1.1 (Git: 10/02/2015) 


OS: Ubuntu x64 
Python 2.7.6 

Time: 
Epoch timestamp: 1420086600 
Timestamp in milliseconds: 1420086600000 
Human time (GMT): Thu, 01 Jan 2015 04:30:00 GMT 

ISS (ZARYA)    
1 25544U 98067A 15096.52834639 .00016216 00000-0 24016-3 0 9993 
2 25544 51.6469 82.0200 0006014 185.1879 274.8446 15.55408008936880 

observation point: N0 E0 alt=0 

Test 1: 

Gpredict: (Time, Az, El, Slant Range, Range Velocity) 

2015 01 01 04:30:00 202.31 -21.46 5638 -5.646 
2015 01 01 04:40:00 157.31 -2.35 2618 -3.107 
2015 01 01 04:50:00 72.68 -10.26 3731 5.262 

Pyephem 3.7.5.3 (default atmospheric refraction) 

(2015/1/1 04:30:00, 202:18:45.3, -21:27:43.0, 5638.0685, -5.3014228515625) 
(2015/1/1 04:40:00, 157:19:08.3, -1:21:28.6, 2617.9915, -2.934402099609375) 
(2015/1/1 04:50:00, 72:40:59.9, -10:15:15.1, 3730.78375, 4.92381201171875) 
No atmospheric refraction 
(2015/1/1 04:30:00, 202:18:45.3, -21:27:43.0, 5638.0685, -5.3014228515625) 
(2015/1/1 04:40:00, 157:19:08.3, -1:21:28.6, 2617.9915, -2.934402099609375) 
(2015/1/1 04:50:00, 72:40:59.9, -10:15:15.1, 3730.78375, 4.92381201171875) 

Pypredict 

1420086600.0 
{'decayed': 0, 'elevation': -19.608647085869123, 'name': 'ISS (ZARYA)', 'norad_id': 25544, 'altitude': 426.45804846615556, 'orbit': 92208, 'longitude': 335.2203454719759, 'sunlit': 1, 'geostationary': 0, 'footprint': 4540.173580837984, 'epoch': 1420086600.0, 'doppler': 1635.3621339278857, 'visibility': 'D', 'azimuth': 194.02436209048014, 'latitude': -45.784314563471646, 'orbital_model': 'SGP4', 'orbital_phase': 73.46488929141783, 'eclipse_depth': -8.890253049060693, 'slant_range': 5311.3721164183535, 'has_aos': 1, 'orbital_velocity': 27556.552465256085} 
1420087200.0 
{'decayed': 0, 'elevation': -6.757496200551716, 'name': 'ISS (ZARYA)', 'norad_id': 25544, 'altitude': 419.11153234752874, 'orbit': 92208, 'longitude': 9.137628905963876, 'sunlit': 1, 'geostationary': 0, 'footprint': 4502.939901708917, 'epoch': 1420087200.0, 'doppler': 270.6901377419433, 'visibility': 'D', 'azimuth': 139.21315598291235, 'latitude': -20.925997669236732, 'orbital_model': 'SGP4', 'orbital_phase': 101.06301876416072, 'eclipse_depth': -18.410968838249545, 'slant_range': 3209.8444916123644, 'has_aos': 1, 'orbital_velocity': 27568.150821416708} 
1420087800.0 
{'decayed': 0, 'elevation': -16.546383900323555, 'name': 'ISS (ZARYA)', 'norad_id': 25544, 'altitude': 414.1342802649042, 'orbit': 92208, 'longitude': 31.52356804788407, 'sunlit': 1, 'geostationary': 0, 'footprint': 4477.499436144489, 'epoch': 1420087800.0000002, 'doppler': -1597.032808834609, 'visibility': 'D', 'azimuth': 76.1840387294104, 'latitude': 9.316828913183791, 'orbital_model': 'SGP4', 'orbital_phase': 128.66115193399546, 'eclipse_depth': -28.67721196244149, 'slant_range': 4773.838774518728, 'has_aos': 1, 'orbital_velocity': 27583.591664378775} 


Test 2 (short time): 
Gpredict: (Slant Range, Range Velocity) 
2015 01 01 04:30:00 5638 -5.646 
2015 01 01 04:30:10 5581 -5.648 
->5.7 km/s avg 

(2015/1/1 04:30:00, 5638.0685, -5.3014228515625) 
(2015/1/1 04:30:10, 5581.596, -5.30395361328125) 
->5.7 km/s avg 

Pyephem

import ephem 
import time 
#TLE Kepler elements 

line1 = "ISS (ZARYA)"    
line2 = "1 25544U 98067A 15096.52834639 .00016216 00000-0 24016-3 0 9993" 
line3 = "2 25544 51.6469 82.0200 0006014 185.1879 274.8446 15.55408008936880" 

satellite = ephem.readtle(line1, line2, line3) # create ephem object from tle information 

obs = ephem.Observer() # recreate Oberserver with current time 
obs.lon, obs.lat, obs.elevation = '0' , '0' , 0 

print('Pyephem Default (atmospheric refraction)') 

obs.date = '2015/1/1 04:30:00' 
satellite.compute(obs) 
print(obs.date, satellite.az, satellite.alt,satellite.range/1000, satellite.range_velocity/1000) 

obs.date = '2015/1/1 04:40:00' 
satellite.compute(obs) 
print(obs.date, satellite.az, satellite.alt,satellite.range/1000, satellite.range_velocity/1000) 

obs.date = '2015/1/1 04:50:00' 
satellite.compute(obs) 
print(obs.date, satellite.az, satellite.alt,satellite.range/1000, satellite.range_velocity/1000) 

obs.pressure = 0 # disable atmospheric refraction 
print('Pyephem No atmospheric refraction') 

obs.date = '2015/1/1 04:30:00' 
satellite.compute(obs) 
print(obs.date, satellite.az, satellite.alt,satellite.range/1000, satellite.range_velocity/1000) 

obs.date = '2015/1/1 04:40:00' 
satellite.compute(obs) 
print(obs.date, satellite.az, satellite.alt,satellite.range/1000, satellite.range_velocity/1000) 

obs.date = '2015/1/1 04:50:00' 
satellite.compute(obs) 
print(obs.date, satellite.az, satellite.alt,satellite.range/1000, satellite.range_velocity/1000) 

print('10 s timing') 
obs.date = '2015/1/1 04:30:00' 
satellite.compute(obs) 
print(obs.date, satellite.range/1000, satellite.range_velocity/1000) 
obs.date = '2015/1/1 04:30:10' 
satellite.compute(obs) 
print(obs.date, satellite.range/1000, satellite.range_velocity/1000) 

Pypredict

import predict 
import datetime 
import time 

format = '%Y/%m/%d %H:%M:%S' 

tle = """ISS (ZARYA) 
1 25544U 98067A 15096.52834639 .00016216 00000-0 24016-3 0 9993 
2 25544 51.6469 82.0200 0006014 185.1879 274.8446 15.55408008936880""" 
qth = (0, 10, 0) # lat (N), long (W), alt (meters) 

#expect time as epoch time float 

time= (datetime.datetime.strptime('2015/1/1 04:30:00', format) -datetime.datetime(1970,1,1)).total_seconds() 
result = predict.observe(tle, qth, time) 
print time 
print result 

time= (datetime.datetime.strptime('2015/1/1 04:40:00', format) -datetime.datetime(1970,1,1)).total_seconds() 
result = predict.observe(tle, qth, time) 
print time 
print result 

time= (datetime.datetime.strptime('2015/1/1 04:50:00', format) -datetime.datetime(1970,1,1)).total_seconds() 
result = predict.observe(tle, qth, time) 
print time 
print result 

отладки выход Gpredict и PyEphem

PyPredict

Name = ISS (ZARYA) 
current jd = 2457023.68750 
current mjd = 42003.7 
satellite jd = 2457119.02835 
satellite mjd = 42099 
SiteLat = 0 
SiteLong = 6.28319 
SiteAltitude = 0 
se_EPOCH : 115096.52834638999775052071 
se_XNO :   0.06786747737871574870 
se_XINCL :   0.90140843391418457031 
se_XNODEO :   1.43151903152465820312 
se_EO  :   0.00060139998095110059 
se_OMEGAO :   3.23213863372802734375 
se_XMO :   4.79694318771362304688 
se_BSTAR :   0.00024016000679694116 
se_XNDT20 :   0.00000000049135865048 
se_orbit :       93688 
dt  : -137290.81880159676074981689 
CrntTime = 42004.2 
SatX = -3807.5 
SatY = 2844.85 
SatZ = -4854.26 
Radius = 6793.68 
SatVX = -5.72752 
SatVY = -3.69533 
SatVZ = 2.32194 
SiteX = -6239.11 
SiteY = 1324.55 
SiteZ = 0 
SiteVX = -0.0965879 
SiteVY = -0.454963 
Height = 426.426 
SSPLat = -0.795946 
SSPLong = 0.432494 
Azimuth = 3.53102 
Elevation = -0.374582 
Range = 5638.07 
RangeRate = -5.30142 
(2015/1/1 04:30:00, 5638.0685, -5.3014228515625) 

Gpredict 

time: 2457023,687500 
pos obs: -6239,093574, 1324,506494, 0,000000 
pos sat: -3807,793748, 2844,641722, -4854,112635 
vel obs: -0,096585, -0,454962, 0,000000 
vel sat: -6,088242, -3,928388, 2,468585 

Gpredict (sgp_math.h)

/------------------------------------------ ------------------------/

/* Converts the satellite's position and velocity */ 
/* vectors from normalised values to km and km/sec */ 
void 
Convert_Sat_State(vector_t *pos, vector_t *vel) 
{ 
     Scale_Vector(xkmper, pos); 
     Scale_Vector(xkmper*xmnpda/secday, vel); 

} /* Procedure Convert_Sat_State */ 

Ephem (Libastro)

*SatX = ERAD*posvec.x/1000; /* earth radii to km */ 
*SatY = ERAD*posvec.y/1000; 
*SatZ = ERAD*posvec.z/1000; 
*SatVX = 100*velvec.x;  /* ?? */ 
*SatVY = 100*velvec.y; 
*SatVZ = 100*velvec.z; 
0

Обновление до самого последнего выпуска из pyephem (я пробовал V3.7.6.0), похоже, решает проблему. Диапазон скорости теперь точно соответствует значениям, данным другим обычно используемым программным обеспечением для отслеживания.