利用python通過GPS的N文件計算所有衛星的軌道坐標

02-23

接上一篇

原理與GPS的RINEX文件格式百度既可以知道，希望我的代碼會給剛學開發的學弟學妹們點啟發，如有不足之處也歡迎指正，共同交流！

直接代碼如下：

# -*- coding: utf-8 -*-"""Created on Tue Dec 12 21:09:54 2017@author: Administrator"""import mathfrom numpy import *import calendarc_0=299792458 #光速with open(rC:UsersAdministratorDesktopGPS試驗田DATA試驗數據d1jz_O.txt,r,encoding=utf-8) as f1: ofile_lines=f1.readlines() #按行讀取O文件with open(rC:UsersAdministratorDesktopGPS試驗田DATA試驗數據d1jz_N.txt,r,encoding=utf-8) as f2: nfile_lines=f2.readlines() #按行讀取N文件#print(x[1])#with open(rC:UsersAdministratorDesktopGPS試驗田DATA軌道坐標.txt,w,encoding=utf-8) as f3:# f3.write() #初始化軌道坐標文檔########對O文件中所有曆元厲遍求其對應衛星軌道坐標#######for ly_a in ofile_lines: odic_in_file=eval(ly_a)# print(odic_in_file) liyuan_o=odic_in_file[曆元] satellite_nums=[] for i in range(odic_in_file[衛星數量]): sl_num=odic_in_file[衛星+str(i+1)] satellite_nums.append(sl_num) #print(satellite_nums) #print(liyuan) liyuan_sl=liyuan_o.split() #print(liyuan_sl) #week_sc=float(liyuan_sl[] week_num=calendar.weekday(int(20+liyuan_sl[0]),int(liyuan_sl[1]),int(liyuan_sl[2]))+1 #調用python日曆應用，直接讀取星期，但python里周一的代碼是0,所以要進行重新設計下 if week_num==7: week_num=0 #print(week_num) #t_c_dl= week_sc=(week_num*24+int(liyuan_sl[3]))*3600+int(liyuan_sl[4])*60+float(liyuan_sl[5]) #求出周內秒 (帶有傳輸時間延時的) #####求每顆衛星的發射信號時間（去掉接收機接收時間延時）####### satellites_tx_list=[] for esa in satellite_nums: t_c1=odic_in_file[esa+號衛星C1] t_p2=odic_in_file[esa+號衛星P2] tx1=t_c1/c_0 tx2=t_p2/c_0 week_sc_c1=week_sc-tx1 week_sc_p2=week_sc-tx2 tx_list=[week_sc_c1,week_sc_p2] satellites_tx_list.append(tx_list) #最終衛星軌道的發射信號時的曆元的周內秒，列表中的列表，裡面列表裡0是c1,1是p2 #print(satellites_tx_list) #print(week_sc-439200.0) ########對比曆元####### ndic_in_file=eval(nfile_lines[0]) liyuan_n=ndic_in_file[曆元] liyuan_sln=liyuan_n.split() try: int(20+liyuan_sl[0])==int(20+liyuan_sln[0]) and int(liyuan_sl[1])==int(liyuan_sln[1]) and int(liyuan_sl[2])==int(liyuan_sln[2])# print(對上了) #判斷其曆元是否在同一天內 except: print(O文件與N文件曆元不符) #print(liyuan_sln) ######獲取最小曆元差值的組數###### stl_num=int(satellite_nums[0]) dai_min_dic_list=[] for dic_n_s in nfile_lines: dic_nn=eval(dic_n_s) if dic_nn[衛星PRN號]==stl_num: dai_min_dic_list.append(dic_nn) #print(dai_min_dic_list) #print(type(dai_min_dic_list[0])) dai_min_list=[] for ii in dai_min_dic_list: dai_liyuan_x=ii[曆元] dai_liyuan=dai_liyuan_x.split() shifen=abs(int(dai_liyuan[3])*60+int(dai_liyuan[4])-int(liyuan_sl[3])*60-int(liyuan_sl[4])) dai_min_list.append(shifen) min_num_position=dai_min_list.index(min(dai_min_list)) #最小差值的位置 #print(min_num_position) ####################################################################################### ############## 計算O文件對應曆元的衛星軌道坐標 ####################################### ####################################################################################### def satellite_coord(t,n,sqrt_A,TEO,M0,e,w,C_us,C_uc,C_rs,C_rc,C_is,C_ic,I_0,i,OMEGA,OMEGA_DOT): GM=3.986005e14 a=math.pow(sqrt_A,2) w_ie=7.292115e-5 #######以下為未知數######## n_x=n t_oe=TEO M_0=M0 e=e w=w C_us=C_us C_uc=C_uc C_rs=C_rs C_rc=C_rc C_is=C_is C_ic=C_ic i_0=I_0 i_top_2=i omega_0=OMEGA omega_dian=OMEGA_DOT ############基本公式################## n_0=math.sqrt(GM/math.pow(a,3)) #計算理論角速度 a為軌道長軸，N文件給出 n=n_0+n_x # n為平均角速度 n_0 為常數理論角速度 n_x為平均運行速度差，N文件提供 t_x=t-t_oe ####################t_x為歸化時間 t為接收機收到該電文時接收機的GPS時間（是整個公式中唯一唯一唯一的變數！！！！！！！） t_oe為參考考曆元，由N文件提供 M=M_0+n*t_x #計算 M為曆元t的平近角點(單位弧度) m_0的值在N文件中給出 n 和t_x 上文計算出 #E=M+e*math.sin(E) ##########需要用迭代法來求解！這個待定！！！！比較難理解！！需要單獨定義個迭代函數！！#########E為偏近點角(單位弧度) e為軌道偏心率 N文件中給出 E_dic={} ##將迭代放入一個字典中進行，便於理解和命名迭代法算E值 E_dic[E0]=M for i in range(100000): E_dic[E+str(i+1)]=M+e*math.sin(E_dic[E+str(i)]) if abs(E_dic[E+str(i+1)]-E_dic[E+str(i)])<0.0000000001: E=E_dic[E+str(i+1)]# print(i)# print(E2值為：+str(E)) break # print(E值為：+str(E)) r_0=a*(1-e*math.cos(E)) #r_0為衛星的地心矢徑 # print(r_0) # print(math.sqrt((1+e)/(1-e))*math.tan(E/2):+str((1+e)/(1-2)*math.tan(E/2))+math.tan(E/2):+str(math.tan(E/2))) f=2*math.atan(math.sqrt((1+e)/(1-e))*math.tan(E/2)) #f為真近點角 q_0=w+f #q_0 為升角距，w為軌道近地點角距，由N文件提供 b_u=C_us*math.sin(2*q_0)+C_uc*math.cos(2*q_0) #b_u b_r b_i 為升交點角距、地心失徑、軌道面傾角攝動改正項 C_us......C_ic 由N文件提供 b_r=C_rs*math.sin(2*q_0)+C_rc*math.cos(2*q_0) b_i=C_is*math.sin(2*q_0)+C_ic*math.cos(2*q_0) q=q_0+b_u #q為經過改正後的升交點角距 r_shi=r_0+b_r #r_shi為矢量，是衛星失徑 i=i_0+b_i+i_top_2*t_x #i為衛星軌道傾角 i_0為N文件提供 i_top_2標識i的兩階導數，是表示加速度的意思，這裡的值有N文件提供，t_x上文給出 kai=omega_0+(omega_dian-w_ie)*t_x-w_ie*t_oe #kai 表示希臘字母「入」，這裡的意義是曆元t的升交點的經度，omega_0為該時刻的升交點赤經，omega_dian為赤經變化率,這兩個都為N文件提供，w_ie為地球角速率，已知量，固定值，t_oe為參考時刻的曆元，已知量，N文件提供 #####衛星在軌道直角坐標系中的坐標############ coord_0_mat=r_shi*mat([[math.cos(q)],[math.sin(q)],[0]]) #x_0 y_0 z_0是衛星在軌道直角坐標系中的坐標 #########衛星在協議地球坐標系的直角坐標###### coord_t_mat=mat([[math.cos(kai),-math.sin(kai)*math.cos(i),math.sin(kai)*math.sin(i)],[math.sin(kai),math.cos(kai)*math.cos(i),-math.cos(kai)*math.sin(i)],[0,math.sin(i),math.cos(i)]])*coord_0_mat ###########坐標輸出########### # print(協議坐標系中坐標為： +X=+str(coord_t_mat[0,0])+ +Y=+str(coord_t_mat[1,0])+ +Z=+str(coord_t_mat[2,0])) coord_list=[coord_t_mat[0,0],coord_t_mat[1,0],coord_t_mat[2,0]] return coord_list dic_zuobiao={} #現在一個小循環里測試，最後還是以列表形式，寫入文件進行存檔 count_num=0 for sa in range(len(satellite_nums)): t_c1_li=satellites_tx_list[count_num][0] t_p2_li=satellites_tx_list[count_num][1] count_num+=1 stl_num=int(satellite_nums[sa]) min_dic_list=[] for dic_n_s in nfile_lines: dic_nn=eval(dic_n_s) if dic_nn[衛星PRN號]==stl_num: min_dic_list.append(dic_nn)# print(min_num_position)# print(min_dic_list) #容易出現各種list超限的情況，所以要用來調試 try: daisuan_dic=min_dic_list[min_num_position] #最終要進行計算的衛星所對應的曆元的數據（字典形式） (PS:由於有特殊變態情況就是某顆衛星的狀態可能就又一個曆元的軌道參數或者曆元不足，這跟觀測條件有關，所以需要異常的判斷) except: for i in range(min_num_position): gezhongshi=min_num_position-1-i try: daisuan_dic=min_dic_list[gezhongshi] break except: continue # print(daisuan_dic[sqrt_A]) dic_zuobiao[str(stl_num)+號衛星+liyuan_o+C1碼求得坐標]=satellite_coord(t_c1_li,daisuan_dic[n],daisuan_dic[sqrt_A],daisuan_dic[TEO],daisuan_dic[M0],daisuan_dic[e],daisuan_dic[w],daisuan_dic[C_us],daisuan_dic[C_uc],daisuan_dic[C_rs],daisuan_dic[C_rc],daisuan_dic[C_is],daisuan_dic[C_ic],daisuan_dic[I_0],daisuan_dic[i],daisuan_dic[OMEGA],daisuan_dic[OMEGA_DOT]) dic_zuobiao[str(stl_num)+號衛星+liyuan_o+P2碼求得坐標]=satellite_coord(t_p2_li,daisuan_dic[n],daisuan_dic[sqrt_A],daisuan_dic[TEO],daisuan_dic[M0],daisuan_dic[e],daisuan_dic[w],daisuan_dic[C_us],daisuan_dic[C_uc],daisuan_dic[C_rs],daisuan_dic[C_rc],daisuan_dic[C_is],daisuan_dic[C_ic],daisuan_dic[I_0],daisuan_dic[i],daisuan_dic[OMEGA],daisuan_dic[OMEGA_DOT])# print(dic_zuobiao)# print(dic_zuobiao) with open(rC:UsersAdministratorDesktopGPS試驗田DATA軌道坐標.txt,a,encoding=utf-8) as f4: f4.write(str(dic_zuobiao)+ )