本文主要是介绍根据广播星历计算GNSS卫星在瞬时地球坐标系中的坐标，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！
由广播星历计算GNSS卫星坐标
- 前言
- 读取文件
- 计算卫星坐标
- 可视化
- 静态文本变成透明的（通过类向导添加）：
- 执行如下操作改变背景
- 小结
- 补充
前言
本文是一个简单的由广播星历计算GNSS卫星坐标的程序，可以作为一个仅用于计算坐标的程序使用，先将初始版码在这里，之后会对该程序进行一些优化。
附：这篇文章主要是为了梳理自己的一些想法、思路，高手忽略~不过里面有一些细节问题可能可以给大家带来一些帮助。
注：但是发现了一个很麻烦的情况，就是不管怎么算，都会出现几个数超过了10m的误差，之后好好回想了一下，个人感觉可能是时间参数出了一点点问题，因为程序中没有用到时间改正模型，但是由于GNSS结课后本人就把相关文件都清空了，调试起来比较麻烦，感兴趣的朋友自己试试吧！
读取文件
计算卫星坐标的第一步是获取广播星历中的数据。本文以GPS导航电文文件为例读取文件。
如图，在文件头之后的就是我们需要的数据，数据结构是每个星历块八行，第一行十个数据，之后每行四个。数据的含义比较复杂，可以参考其他文献。网上应该很多相关资料。
读取文件后，数据储存在一个结构体指针中。
struct EPHEMERISBLOCK//每小时一个卫星对应一个基本星历块
{//PRN号 int PRN;double a0, a1, a2;//时间改正数//六个轨道参数double IODE, Crs, Deltan, M0;// ORBIT - 1double Cuc, e, Cus, SqrtA;// ORBIT - 2double Toe, Cic, OMEGA, Cis;// ORBIT - 3double i0, Crc, omega, OMEGAdot;// ORBIT - 4double IDOT, GpsWeekNumber, L2C, L2P;// ORBIT - 5double SatAccuracy, SatHealth, TGD, IODC;// ORBIT - 6
};##
构造一个计算类CMyCalculation，专门用于完成读取文件和计算功能,在类中添加以下变量和方法：
private:int EphemerisBlockNum;//记录星历块数目EPHEMERISBLOCK *m_pGpsEphemeris;//记录星历数据public void ReadBrodcastEphemeris(CString strEpheNam);
以下为读取文件函数的实现：
void CMyCaculation::ReadBrodcastEphemeris(CString strEpheNam)
{int HeadLineNum = 0;int WeekNo;double WeekSecond;//打开文件CStdioFile pfEph;BOOL IsEn = pfEph.Open(strEpheNam, CFile::modeRead);if (!IsEn) {AfxMessageBox(_T("文件读取失败"),MB_OK,0);return;}//读入头文件CString strLine;while (IsEn){IsEn = pfEph.ReadString(strLine);HeadLineNum++;int index = strLine.Find(CString("END OF HEADER"));if (-1 != index)break;}//计算星历块数int AllNum = 0;while (IsEn){IsEn = pfEph.ReadString(strLine);AllNum++;}//指针不为空时释放指针内存，之后重新赋值if(NULL!=m_pGpsEphemeris){delete m_pGpsEphemeris;}//临时读入星历块EphemerisBlockNum = (AllNum + 1) / 8;m_pGpsEphemeris = new EPHEMERISBLOCK[EphemerisBlockNum];GPSTIME *pGpsTime = new GPSTIME[EphemerisBlockNum];if (!m_pGpsEphemeris || !pGpsTime) return;//将文件指针调整到数据位置pfEph.SeekToBegin();for (int i = 0; i<HeadLineNum; i++)IsEn = pfEph.ReadString(strLine);//定义读取的参数int mPrn;//卫星号PRNo int year, month, day, hour, minute;//卫星钟参考时刻double msecond;double a0, a1, a2;//卫星钟飘参数double IODE, Crs, DeltN, M0;//数据星历发布时间，在轨道径向方向上周期改正正弦的振幅double Cuc, e, Cus, sqrtA;//轨道延迹方向上周期改正余弦振幅 、扁心率、轨道延迹方向上周期改正正弦振幅 、长半轴平方根 double Toe, Cic, OMEGA, Cis;//星历参考时刻、轨道倾角周期改正余弦项振幅、参考时刻升交点赤径主项、轨道倾角周期改正正弦项振幅double i0, Crc, omega, OMEGADOT;//参考时间轨道倾角、在轨道径向方向上周期改正余余弦的振幅、近地点角距、升交点赤径在赤道平面中的长期变化double IDOT, L2C, GPSWeek, L2P;轨道倾角变化率、？？、gps周double AccuracyofSat, HealthofSat, TGD, IODC;//卫星精度、卫星健康、电离层群迟改正数//将以下循环中读取数据的CString都改为char*类型strlinefor (int i = 0; i<EphemerisBlockNum; i++){//读取卫星PRN号，星历参考时间IsEn = pfEph.ReadString(strLine);strLine.Replace('D', 'e');char* strline = CSTR2CAHR(strLine);//sscanf_s(strline, "%d %d %d %d %d %d %lf %lf %lf %lf", &mPrn, &year, &month, &day, &hour, &minute, &msecond, &a0, &a1, &a2);year += 2000;WeekNo = Calendar2GpsTime(year, month, day, hour, minute, msecond, WeekSecond);pGpsTime[i].weekno = WeekNo;pGpsTime[i].weekSecond = WeekSecond;IsEn = pfEph.ReadString(strLine);strLine.Replace('D', 'e');strline = CSTR2CAHR(strLine);sscanf_s(strline, "%lf %lf %lf %lf", &IODE, &Crs, &DeltN, &M0);//读 Cuc,e,Cus,sqrtAIsEn = pfEph.ReadString(strLine);strLine.Replace('D', 'e');strline = CSTR2CAHR(strLine);sscanf_s(strline, "%lf %lf %lf %lf", &Cuc, &e, &Cus, &sqrtA);//Toe,Cic,OMEGA,Cis;IsEn = pfEph.ReadString(strLine);strLine.Replace('D', 'e');strline = CSTR2CAHR(strLine);sscanf_s(strline, "%lf %lf %lf %lf", &Toe, &Cic, &OMEGA, &Cis);//i0,Crc,w,OMEGADOTIsEn = pfEph.ReadString(strLine);strLine.Replace('D', 'e');strline = CSTR2CAHR(strLine);sscanf_s(strline, "%lf %lf %lf %lf", &i0, &Crc, &omega, &OMEGADOT);//IDOT,L2Cod,GPSWeek,L2PCodIsEn = pfEph.ReadString(strLine);strLine.Replace('D', 'e');strline = CSTR2CAHR(strLine);sscanf_s(strline, "%lf %lf %lf %lf", &IDOT, &L2C, &GPSWeek, &L2P);//AccuracyofSat,HealthofSat,TGD,IODCIsEn = pfEph.ReadString(strLine);strLine.Replace('D', 'e');strline = CSTR2CAHR(strLine);sscanf_s(strline, "%lf %lf %lf %lf", &AccuracyofSat, &HealthofSat, &TGD, &IODC);//IsEn = pfEph.ReadString(strLine);//赋值m_pGpsEphemeris[i].PRN = mPrn;m_pGpsEphemeris[i].a0 = a0;m_pGpsEphemeris[i].a1 = a1;m_pGpsEphemeris[i].a2 = a2;//&IODE, &Crs, &DeltN, &M0m_pGpsEphemeris[i].IODE = IODE;m_pGpsEphemeris[i].Crs = Crs;m_pGpsEphemeris[i].Deltan = DeltN;m_pGpsEphemeris[i].M0 = M0;//&Cuc, &e, &Cus, &sqrtAm_pGpsEphemeris[i].Cuc = Cuc;m_pGpsEphemeris[i].e = e;m_pGpsEphemeris[i].Cus = Cus;m_pGpsEphemeris[i].SqrtA = sqrtA;//Toe,Cic,OMEGA,Cis;m_pGpsEphemeris[i].Toe = Toe;m_pGpsEphemeris[i].Cic = Cic;m_pGpsEphemeris[i].OMEGA = OMEGA;m_pGpsEphemeris[i].Cis = Cis;//i0,Crc,omega,OMEGADOTm_pGpsEphemeris[i].i0 = i0;m_pGpsEphemeris[i].Crc = Crc;m_pGpsEphemeris[i].omega = omega;m_pGpsEphemeris[i].OMEGAdot = OMEGADOT;//iDOT,L2Cod,GPSWeek,L2PCodm_pGpsEphemeris[i].IDOT = IDOT;m_pGpsEphemeris[i].L2C = L2C;m_pGpsEphemeris[i].L2P = L2P;m_pGpsEphemeris[i].GpsWeekNumber = GPSWeek;//AccuracyofSat,HealthofSat,TGD,IODCm_pGpsEphemeris[i].SatAccuracy = AccuracyofSat;m_pGpsEphemeris[i].SatHealth = HealthofSat;m_pGpsEphemeris[i].TGD = TGD;m_pGpsEphemeris[i].IODC = IODC;}AfxMessageBox(_T("文件读取成功！"), MB_OK, 0);pfEph.Close();if (pGpsTime) { delete[]pGpsTime; pGpsTime = NULL; }}读取文件后，数据存储于结构体中，在需要调用时，只要匹配结构体中的卫星号和用户输入的卫星号就可以获取相应的星历块了。
上述代码中有一个CSTR2CAHR(strLine)方法，这是因为在调试时，笔者发现在使用CString时，出现了CString和const char*不能互相转化的问题，这也是MFC常见的问题，于是定义了一个方法来完成转化：
char* CMyCaculation::CSTR2CAHR(CString str) {const size_t strSize = (str.GetLength() + 1) * 2;char* p = new char[strSize];size_t sz = 0;wcstombs_s(&sz, p, strSize, str, _TRUNCATE);int n = atoi((const char*)p);return p;
}
最终返回一个指针，用于代替strLine存储数据，如：
IsEn = pfEph.ReadString(strLine);strLine.Replace('D', 'e');char* strline = CSTR2CAHR(strLine);//sscanf_s(strline, "%d %d %d %d %d %d %lf %lf %lf %lf", &mPrn, &year, &month, &day, &hour, &minute, &msecond, &a0, &a1, &a2);year += 2000;
经过这些处理，文件中的数据就被存到了m_pGpsEphemeris中，接下来调用其中的数据计算。
计算卫星坐标
添加结构体，用于记录卫星位置。
struct Position {double X;double Y;double Z;};
加入以下方法计算位置：
Position CMyCaculation::CalculateSP(int year, int month, int date, int hour, int minute, int second, int PNR);该函数返回一个Position，便于在之后的结果显示中调用计算结果。
函数的实现代码：
Position SainEarth;double weekSecond;Calendar2GpsTime(year, month, date, hour, minute, second, weekSecond);//调用函数将日历时间转化为周秒//如果获取时间不合法，则跳出函数if (-1 == weekSecond) {AfxMessageBox(_T("输入的时间不合法！"), MB_OK, 0);SainEarth.X= SainEarth.Y=SainEarth.Z = 0;return SainEarth;}//获取用户输入的卫星序列号在数据结构体中的位置int i = 0;//遍历数据块，寻找卫星号相同时的数据for (; i < EphemerisBlockNum; i++) {if (PNR == m_pGpsEphemeris[i].PRN) {break;}}if (PNR != m_pGpsEphemeris[i].PRN) {AfxMessageBox(_T("文件中不包含该卫星的信息！"), MB_OK, 0);}
上述代码中，Calendar2GpsTime为时间转化函数，较为简单就不贴了。
另外，由于个人的疏忽，把指针名字打太长了，直接调用可能不是很好看，就做了一些额外的工作（实际上并不需要这些步骤，如果要美化代码的话可以考虑更改指针名字）。以下为额外工作的内容（并不重要，可以跳过，但为了程序的完整性还是放上来）：
//为之后的计算看起来简洁，先用简单的名字代替结构体中的数据const double GM = 3.986005E14;const double w = 7.292115147E-5;//地球自转角速度double sqrtA = m_pGpsEphemeris[i].SqrtA;double DeltN = m_pGpsEphemeris[i].Deltan;double M0 = m_pGpsEphemeris[i].M0;double Toe = m_pGpsEphemeris[i].Toe;double e = m_pGpsEphemeris[i].e;double ws = m_pGpsEphemeris[i].omega;double Cuc = m_pGpsEphemeris[i].Cuc;double Crc = m_pGpsEphemeris[i].Crc;double Cic = m_pGpsEphemeris[i].Cic;double Cus = m_pGpsEphemeris[i].Cus;double Crs = m_pGpsEphemeris[i].Crs;double Cis = m_pGpsEphemeris[i].Cis;double di = m_pGpsEphemeris[i].IDOT;double i0 = m_pGpsEphemeris[i].i0;double OMEGA = m_pGpsEphemeris[i].OMEGA;double OMEGADOT = m_pGpsEphemeris[i].OMEGAdot;double t = weekSecond;
这些变量的含义在之后的代码中会有必要的注释说明。
之后定义计算过程中必须使用的变量：
private: double n, Ms, r0, fs, u0;//平均角速度，平近点角，近地点角距，真近点角double deltU, deltR, deltI;//摄动参数double U, R, I;//改正后的参数double x0, y0, z0;//卫星在轨道坐标系中的坐标
考虑到程序持续工作的可能，在类中声明可能是一种比较好的选择。
之后就可以计算坐标了，结果储存在SainEarth中：
n = sqrt(GM / pow(sqrtA, 6)) + DeltN;//GM=3.986005E14,计算平均角速度Ms = M0 + n*(t - Toe);//计算平近点角double Es = 0.0; //偏近点角double preEs = 0.0;do {preEs = Es;Es = Ms + e*sin(preEs);} while (Es - preEs >= 10E-12);//计算偏近点角//计算卫星的近地点角距r0 = pow(sqrtA, 2)*(1 - e*cos(Es));//a=sqrtA^2;//计算真近点角，需要根据象限来调整真近点角//fs = atan(sqrt(1.0 - pow(e,2))*sin(Es) / (cos(Es) - e));double thetax = (cos(Es) - e);double thetay = sqrt(1.0 - pow(e, 2))*sin(Es);if (thetay >= 0 && thetax >= 0) {fs = atan(thetay / thetax);}else if (thetay >= 0 && thetax < 0) {fs = PI +atan(thetay / thetax);}else if (thetay < 0 && thetax < 0) {fs = PI + atan(thetay / thetax);}else {fs = 2 * PI + atan(thetay / thetax);}//计算升交距角u0u0 = ws + fs; //ws为近地点角距，即omega//计算摄动改正项deltU = Cuc*cos(2 * u0) + Cus*sin(2 * u0);deltR = Crc*cos(2 * u0) + Crs*sin(2 * u0);deltI = Cic*cos(2 * u0) + Cis*sin(2 * u0);//改正相关参数U = u0 + deltU;R = r0 + deltR;I = i0 + di*(t - Toe) + deltI;//i0为参考时间角度，di为轨道倾角变化率，即IDOT//计算卫星在轨道面坐标系中的位置x0 = R*cos(U);y0 = R*sin(U);z0 = 0;//计算观测时刻升交点的经度Ldouble L = OMEGA + OMEGADOT*(t - Toe) - w*t;//计算瞬时地球坐标系中卫星的坐标并输出SainEarth.X = x0*cos(L) - y0*cos(I)*sin(L);SainEarth.Y = x0*sin(L) + y0*cos(I)*cos(L);SainEarth.Z = y0*sin(I);return SainEarth;
到这里，卫星位置的计算就完成了。
可视化
先看看最终的样子吧：
总体来说是还可以的。
这是一个基于对话框的MFC应用，创建对话框并没有什么需要注意的。
界面设计如下:
首先添加一些必要的变量：
CMyCaculation caculation;//在一次读取文件操作中更新CString Prn, Year, Month, Date, Hour, Minute, Second;//用于存储用户输入的信息//CStatic* m_static = (CStatic*)GetDlgItem(IDC_STATIC);//需要改IDint StatueOfFile = 0;//记录文件读取状态，如果成功读入文件，则变为1int TableHead = 0;//记录表头初始化状态，如果不为零，则不再初始化添加 导入星历文件 按钮响应事件：
void CVIsualGNSSDlg::OnBnClickedImport()
{// TODO: 在此添加控件通知处理程序代码//选择文件导入星历文件StatueOfFile = 0;//用于记录星历文件导入状态CString strOpenFileType = _T("(*.*)|*.*|"), fname;CFileDialog FileDlg(TRUE, _T("*.*"), NULL, OFN_HIDEREADONLY | OFN_OVERWRITEPROMPT, strOpenFileType);if (FileDlg.DoModal() == IDOK)fname = FileDlg.GetFolderPath() + _T("\\") + FileDlg.GetFileName();caculation.ReadBrodcastEphemeris(fname);StatueOfFile = 1;
}
添加 计算卫星坐标 按钮事件：
void CVIsualGNSSDlg::OnBnClickedCalculation()
{// TODO: 在此添加控件通知处理程序代码if (StatueOfFile == 0) {AfxMessageBox(_T("未导入星历文件"), MB_OK, 0);return;}Position A=caculation.CalculateSP(_ttoi(Year), _ttoi(Month), _ttoi(Date), _ttoi(Hour), _ttoi(Minute), _ttoi(Second), _ttoi(Prn));//计算坐标//CListCtrl *pList = (CListCtrl*)GetDlgItem(IDC_LIST1);CString Result;CString X,Y,Z;//将坐标转化为宽字符,似乎不需要X.Format(_T("%.10lf"), A.X);Y.Format(_T("%.10lf"), A.Y);Z.Format(_T("%.10lf"), A.Z);Result.Format(_T("%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s"), Prn,Year,Month,Date,Hour,Minute,Second,X, Y, Z);//MessageBox(Result, _T("程序运行结果"), MB_OK);//pList->SetWindowTextW(Result);//计算较差CListBox *pList = (CListBox*)GetDlgItem(IDC_LIST1);if (0==TableHead) {//初始化表格pList->AddString(_T("PRN\tYear\tMonth\tDate\tHour\tMinute\tSecond\tX\t\t\tY\t\t\tZ"));TableHead = 1;}pList->AddString(Result);SetHScroll(pList);
}
添加文本编辑框更改事件：
void CVIsualGNSSDlg::OnEnChangeEdits()
{int ID = LOWORD(GetCurrentMessage()->wParam);CEdit* pBox;pBox = (CEdit*)GetDlgItem(ID);switch (ID) {case IDC_PRN:pBox->GetWindowText(Prn);break;case IDC_YEAR:pBox->GetWindowText(Year);break;case IDC_month:pBox->GetWindowText(Month);break;case IDC_date:pBox->GetWindowText(Date);break;case IDC_hour:pBox->GetWindowText(Hour);break;case IDC_minute:pBox->GetWindowText(Minute);break;case IDC_second:pBox->GetWindowText(Second);break; }看了一下这个函数，感觉血亏，文本框每改变一次都要调用一次这个函数，不过因为数据并不多，也没对程序造成多大影响。之后的更新会考虑把用户数据的更新放在按钮事件里，或者说在设置一个更容易监控的量来决定是否更新数据，减少这个函数的调用。
到这里呢，程序基本上是完成了。由于这个程序不支持最大化，加上本来有一步与精密星历的较差计算，所以ListBox显得有点不够看了，就加了一个ScrollBar。
void CVIsualGNSSDlg::SetHScroll(CListBox *pList)
{CDC* dc = GetDC();SIZE s;int index;CString str;long temp;for(index= 0; index< pList->GetCount(); index++){pList->GetText(index,str);s = dc->GetTextExtent(str,str.GetLength()+1); // 获取字符串的像素大小// IDC_LISTBOX为m_List的资源IDtemp = (long)SendDlgItemMessage(IDC_LIST1, LB_GETHORIZONTALEXTENT, 0, 0); //temp得到滚动条的宽度if (s.cx > temp) {SendDlgItemMessage(IDC_LIST1, LB_SETHORIZONTALEXTENT, (WPARAM)s.cx*3, 0);}}ReleaseDC(dc);
}
这一段实现是有点问题的，可能是因为计算结果太长了，而第一行的数据又比较简短，所以调用时乘了个3。
然后就是一些微调啦~
静态文本变成透明的（通过类向导添加）：
HBRUSH CVIsualGNSSDlg::OnCtlColor(CDC* pDC, CWnd* pWnd, UINT nCtlColor)
{HBRUSH hbr = CDialog::OnCtlColor(pDC, pWnd, nCtlColor);/*if (pWnd->GetDlgCtrlID() == IDC_STATIC){pDC->SetTextColor(RGB(255, 0, 0));}*/if (nCtlColor == CTLCOLOR_STATIC) //静态文本{pDC->SetBkMode(TRANSPARENT); //设置控件透明return (HBRUSH)::GetStockObject(NULL_BRUSH); //记住一定要有这句}return hbr;
}
执行如下操作改变背景
1.点击添加资源
2.导入BMP
3.选择自己想要的图片
上面这些截图可能有点糊，如果影响了理解可以cue我。
到选图片这个环节我就发现图片文件不显示，不管是不是放在res中。于是只好把文件类型改成了全部。
4.更改OnPaint函数：
void CVIsualGNSSDlg::OnPaint()
{if (IsIconic()){CPaintDC dc(this); // 用于绘制的设备上下文SendMessage(WM_ICONERASEBKGND, reinterpret_cast<WPARAM>(dc.GetSafeHdc()), 0);// 使图标在工作区矩形中居中int cxIcon = GetSystemMetrics(SM_CXICON);int cyIcon = GetSystemMetrics(SM_CYICON);CRect rect;GetClientRect(&rect);int x = (rect.Width() - cxIcon + 1) / 2;int y = (rect.Height() - cyIcon + 1) / 2;// 绘制图标dc.DrawIcon(x, y, m_hIcon);}else{//CDialog::OnPaint(); CPaintDC dc(this);CRect rect;GetClientRect(&rect);CDC dcMem;dcMem.CreateCompatibleDC(&dc);CBitmap bmpBackground;bmpBackground.LoadBitmap(IDB_BITMAP2);BITMAP bitmap;bmpBackground.GetBitmap(&bitmap);CBitmap *pbmpOld = dcMem.SelectObject(&bmpBackground);dc.BitBlt(0, 0, rect.Width(), rect.Height(), &dcMem, 0, 0,SRCCOPY);}
}
这段代码可以参考这篇文章：https://blog.csdn.net/gzywh/article/details/6467468?spm=5176.100239.blogcont26210.4.pK8kD4
到此，这个应用程序就算是比较完善的了。
小结
总的来看，这个程序的考虑比较周全，避开了一些可能的错误使用，如用户先点计算而不是导入等问题。可以优化的地方也很多。
以后有空来把这些坑填掉，不过可能就是以其他的形式了。
补充
有同学问到时间转化函数怎么写，其实这个只要基础语法会应该就没什么问题，之前主要是考虑到篇幅过长久没贴，为了照顾新手还是补上了：
int CMyCaculation::Calendar2GpsTime(int nYear, int nMounth, int nDay, int nHour, int nMinute, double dSecond, double &WeekSecond)
{int DayofMonth = 0;int DayofYear = 0;int weekno = 0;int dayofweek;int m;if (nYear < 1980 || nMounth < 1 || nMounth > 12 || nDay < 1 || nDay > 31 || nHour < 0 || nMinute < 0 || dSecond < 0) {return -1;}//计算从1980年到当前的前一年的天数for (m = 1980; m < nYear; m++){if ((m % 4 == 0 && m % 100 != 0) || (m % 400 == 0)){DayofYear += 366;}elseDayofYear += 365;}//计算当前一年内从元月到当前前一月的天数for (m = 1; m < nMounth; m++){if (m == 1 || m == 3 || m == 5 || m == 7 || m == 8 || m == 10 || m == 12)DayofMonth += 31;else if (m == 4 || m == 6 || m == 9 || m == 11)DayofMonth += 30;else if (m == 2){if ((nYear % 4 == 0 && nYear % 100 != 0) || (nYear % 400 == 0))DayofMonth += 29;elseDayofMonth += 28;}}DayofMonth = DayofMonth + nDay - 6;//加上当月天数/减去1980年元月的6日 weekno = (DayofYear + DayofMonth) / 7;//计算GPS周dayofweek = (DayofYear + DayofMonth) % 7;//计算GPS 周秒时间WeekSecond = dayofweek * 86400 + nHour * 3600 + nMinute * 60 + dSecond;return weekno;
}
这篇关于根据广播星历计算GNSS卫星在瞬时地球坐标系中的坐标的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！
