一、北斗系统概述

1.空间星座

北斗卫星导航系统简称北斗系统，英文缩写为 BDS，其空间星座由 5 颗地球静止轨道（GEO）卫星、27 颗中圆地球轨道（MEO）卫星和 3 颗倾斜地球同步轨道（IGSO）卫星组成。GEO 卫星轨道高度 35786 千米，分别定点于东经 58.75 度、80 度、110.5 度、140 度和 160 度；MEO 卫星轨道高度 21528 千米，轨道倾角 55 度；IGSO 卫星轨道高度 35786 千米，轨道倾角 55 度。

2.坐标系统

北斗系统采用 2000 中国大地坐标系（CGCS2000）。CGCS2000 大地坐标系的定义如下：原点位于地球质心； Z 轴指向国际地球自转服务组织（IERS）定义的参考极（IRP）方向； X 轴为 IERS 定义的参考子午面（IRM）与通过原点且同 Z 轴正交的赤道面的交线； Y 轴与 Z、X 轴构成右手直角坐标系。 CGCS2000 原点也用作 CGCS2000 椭球的几何中心，Z 轴用作该旋转椭球的旋转轴。CGCS2000 参考椭球定义的基本常数为：

长半轴：	a = 6378137.0 m
地球（包含大气层）引力常数：
扁率：	f = 1/298.257222101
地球自转角速度：

3.时间系统

北斗系统的时间基准为北斗时（BDT）。BDT 采用国际单位制（SI）秒为基本单位连续累计，不闰秒，起始历元为 2006 年 1 月 1 日协调世界时（UTC）00 时 00 分 00 秒，采用周和周内秒计数。BDT 通过 UTC（NTSC）与国际 UTC 建立联系，BDT 与 UTC 的偏差保持在 100 纳秒以内（模 1 秒）。 BDT 与 UTC 之间的闰秒信息在导航电文中播报。

二、实验目的

卫星星历是描述卫星轨道运动的一组参数，根据卫星星历可计算卫星的实时位置。通过本次实验，熟悉北斗卫星星历形式、基本参数及意义，并能根据卫星星历计算卫星位置。

三、实验内容

利用给定的卫星星历数据（badata.txt ），参考《北斗卫星导航系统空间信号接口控制文件-公开服务信号（2.1 版）》，基于编程实现卫星在 CGCS2000 大地坐标系中的坐标计算。

四、实验过程

卫星的位置是通过卫星历书数据进行计算的。卫星历书数据包含了卫星的轨道参数、时间信息和其他相关参数。通过这些数据，我们可以推算出卫星在特定时间点的准确位置。

1.首先，代码打开一个输入文件，并使用缓冲流将文件内容读取到一个字符串中。这个输入文件包含了卫星历书数据，它的格式可能是一行一组数据，或者其他形式的文本格式。

2.接下来，代码使用链表和哈希映射等数据结构来存储卫星历书数据。链表用于按顺序存储每组数据，而哈希映射则用于快速查找某个特定的卫星数据。每组数据包含了卫星的编号、时间、轨道根数等信息。

3.然后，代码使用存储的数据进行卫星位置的计算。它通过一系列的数学公式和迭代算法来计算卫星的轨道参数。这些轨道参数包括半长轴、偏心率、轨道倾角、升交点赤经等。这些参数是描述卫星轨道形状和位置的重要信息。

4.根据计算得到的轨道参数，代码进一步计算卫星在三维空间中的位置坐标。这个过程涉及到坐标系的转换、时间的计算和向量运算等。

5.最后，代码将计算得到的卫星位置结果写入一个输出文件。输出文件通常包含卫星的编号和对应的位置坐标信息。这样，用户可以通过读取输出文件来获取卫星在特定时间点的准确位置。

6.除了主要的位置计算部分，代码还包括了一些辅助函数。这些函数用于计算儒略日（Julian Day）等时间相关的计算，以及矩阵乘法等基本数学运算。

源码：

import java.io.BufferedReader;
import java.io.FileReader;
import java.io.FileWriter;
import java.io.IOException;
import java.math.BigDecimal;
import java.util.HashMap;
import java.util.LinkedList;public class Test_GNSS {public static void main(String[] args) throws IOException {double miu = 3.986004418e+14;           //CGCS2000 坐标系下的地球引力常数double OmegaEdot = 7.2921150e-5;        //CGCS2000 坐标系下的地球旋转速率double DayToSecond = 24.0 * 3600.0d;    //天秒数//创建字符输入输对象来访问文件FileReader fileReader = new FileReader("D:\\juechen\\badata.txt");FileWriter fileWriter = new FileWriter("D:\\juechen\\badata_result.txt");//创建输入对象的缓冲流BufferedReader bufferedReader = new BufferedReader(fileReader);//将缓冲流写入字符串String inputData = bufferedReader.readLine();//System.out.println(inputData);//创建字符可变的字符串用来记录计算结果StringBuilder outputDta = new StringBuilder("PRN\tx\ty\tz\n");//创建相应的数据结构用来保存原始数据LinkedList<HashMap<String, String>> linkedList = new LinkedList<>();//将卫星星历数据存储在LinkedList中的计算过程for (int i = 0; i < 13; i++) {inputData = bufferedReader.readLine();//将读入的字符串保存在字符串数组中，当匹配到两个及以上的空格时，分割String[] newLineArray = inputData.split("[ ]{2,}");//将列名和所对应的数据存储在HashMap中HashMap<String, String> newLineHashMap = new HashMap<>();newLineHashMap.put("Epoch", newLineArray[0]);newLineHashMap.put("PRN", newLineArray[1]);newLineHashMap.put("WN", newLineArray[2]);newLineHashMap.put("toe", newLineArray[3]);//对超过16位有效位的数用BigDecimal进行有效存储，提高精度BigDecimal db = new BigDecimal(newLineArray[4]);newLineHashMap.put("SqrtA", db.toPlainString());db = new BigDecimal(newLineArray[5]);newLineHashMap.put("e", db.toPlainString());db = new BigDecimal(newLineArray[6]);newLineHashMap.put("i", db.toPlainString());db = new BigDecimal(newLineArray[7]);newLineHashMap.put("Omega0", db.toPlainString());db = new BigDecimal(newLineArray[8]);newLineHashMap.put("w", db.toPlainString());db = new BigDecimal(newLineArray[9]);newLineHashMap.put("M", db.toPlainString());db = new BigDecimal(newLineArray[10]);newLineHashMap.put("deltaN", db.toPlainString());db = new BigDecimal(newLineArray[11]);newLineHashMap.put("OmegaDot", db.toPlainString());db = new BigDecimal(newLineArray[12]);newLineHashMap.put("idot", db.toPlainString());newLineHashMap.put("Crs", newLineArray[13]);newLineHashMap.put("Crc", newLineArray[14]);db = new BigDecimal(newLineArray[15]);newLineHashMap.put("Cus", db.toPlainString());db = new BigDecimal(newLineArray[16]);newLineHashMap.put("Cuc", db.toPlainString());db = new BigDecimal(newLineArray[17]);newLineHashMap.put("Cis", db.toPlainString());db = new BigDecimal(newLineArray[18]);newLineHashMap.put("Cic", db.toPlainString());newLineHashMap.put("toc", newLineArray[19]);db = new BigDecimal(newLineArray[20]);newLineHashMap.put("a0", db.toPlainString());db = new BigDecimal(newLineArray[21]);newLineHashMap.put("a1", db.toPlainString());db = new BigDecimal(newLineArray[22]);newLineHashMap.put("a3", db.toPlainString());//保存原始数据linkedList.add(newLineHashMap);}//用linkList存储的每一组卫星星历数据来计算卫星位置for (int i = 0; i < 13; i++) {HashMap<String, String> calculateHashMap = linkedList.get(i);//计算观测时间的儒略日String[] date_time = calculateHashMap.get("Epoch").split(" ");String date = date_time[0];String time = date_time[1];double[] timeD = new double[]{Double.parseDouble(date.split("-")[0]), Double.parseDouble(date.split("-")[1]), Double.parseDouble(date.split("-")[2])};double[] timeH = new double[]{Double.parseDouble(time.split(":")[0]), Double.parseDouble(time.split(":")[1]), Double.parseDouble(time.split(":")[2])};double EpochT = calculateJD(timeD[0], timeD[1], timeD[2], timeH[0], timeH[1], timeH[2]) - calculateJD(2006d, 1d, 1d, 0, 0, 0);//计算半长轴double A = Math.pow(Double.parseDouble(calculateHashMap.get("SqrtA")), 2);//计算卫星平均角速度double n0 = Math.sqrt(miu / Math.pow(A, 3));//计算观测历元到参考历元的时间差double tk = (EpochT - Double.parseDouble(calculateHashMap.get("WN")) * 7.0d) * DayToSecond - Double.parseDouble(calculateHashMap.get("toe"));//计算改正平均角速度double n = n0 + Double.parseDouble(calculateHashMap.get("deltaN"));//计算平近点角double Mk = Double.parseDouble(calculateHashMap.get("M")) + n * tk;//迭代计算偏近点角double delta = 1.0D;double Ek = Mk;double k = 0L;double e = Double.parseDouble(calculateHashMap.get("e"));while (Math.abs(delta) > 0.0000001D && k <= 1000l) {Ek = Ek - delta / (1 - e * Math.cos(Ek));delta = Ek - e * Math.sin(Ek) - Mk;k = k + 1;}//计算真近点角double sinvk = (Math.sqrt(1 - Math.pow(e, 2)) * Math.sin(Ek)) / (1 - e * Math.cos(Ek));double cosvk = (Math.cos(Ek) - e) / (1 - e * Math.cos(Ek));//计算vkdouble vk = 0;if (sinvk >= 0 && cosvk >= 0) {vk = Math.asin(sinvk);} else if (sinvk >= 0 && cosvk < 0) {vk = Math.PI - Math.asin(sinvk);} else if (sinvk < 0 && cosvk >= 0) {vk = Math.PI * 2 + Math.asin(sinvk);} else if (sinvk < 0 && cosvk < 0) {vk = Math.PI - Math.asin(sinvk);}//计算纬度幅角参数faikdouble faik = vk + Double.parseDouble(calculateHashMap.get("w"));//计算纬度幅角改正项duk、径向改正项drk、轨道倾角改正项dikdouble duk = Double.parseDouble(calculateHashMap.get("Cus")) * Math.sin(2 * faik) + Double.parseDouble(calculateHashMap.get("Cus")) * Math.cos(2 * faik);double drk = Double.parseDouble(calculateHashMap.get("Crs")) * Math.sin(2 * faik) + Double.parseDouble(calculateHashMap.get("Crs")) * Math.cos(2 * faik);double dik = Double.parseDouble(calculateHashMap.get("Cis")) * Math.sin(2 * faik) + Double.parseDouble(calculateHashMap.get("Cis")) * Math.cos(2 * faik);//计算改正后的纬度幅角ukdouble uk = faik + duk;//计算改正后的径向rkdouble rk = A * (1 - e * Math.cos(Ek)) + drk;//计算改正后的轨道倾角ikdouble ik = Double.parseDouble(calculateHashMap.get("i")) + Double.parseDouble(calculateHashMap.get("idot")) * tk + dik;//计算卫星在轨道平面内的坐标(xk,yk)double xk = rk * Math.cos(uk);double yk = rk * Math.sin(uk);//计算历元升交点赤经（地固系）//计算 MEO/IGSO 卫星在 CGCS2000 坐标系中的坐标(Xk,Yk,Zk)double Omegak = 0;double Xk = 0;double Yk = 0;double Zk = 0;if (Double.parseDouble(calculateHashMap.get("i")) > 0) {Omegak = Double.parseDouble(calculateHashMap.get("Omega0")) + (Double.parseDouble(calculateHashMap.get("OmegaDot")) - OmegaEdot) * tk - OmegaEdot * Double.parseDouble(calculateHashMap.get("toe"));Xk = xk * Math.cos(Omegak) - yk * Math.cos(ik) * Math.sin(Omegak);Yk = xk * Math.sin(Omegak) + yk * Math.cos(ik) * Math.cos(Omegak);Zk = yk * Math.sin(ik);}//计算历元升交点的赤经（惯性系）//计算 GEO 卫星在自定义坐标系系中的坐标(Xgk,Ygk,Zgk)//计算 GEO 卫星在 CGCS2000 坐标系中的坐标(Xk,Yk,Zk)double Xgk = 0;double Ygk = 0;double Zgk = 0;double Rxf = 0;double Rzf = 0;double[][] Rz = null;double[][] Rx = null;double[][] Rgk = null;double[][] XYZ = null;if (Double.parseDouble(calculateHashMap.get("i")) == 0) {Omegak = Double.parseDouble(calculateHashMap.get("Omega0")) + Double.parseDouble(calculateHashMap.get("OmegaDot")) * tk - OmegaEdot * Double.parseDouble("toe");Xgk = xk * Math.cos(Omegak) - yk * Math.cos(ik) * Math.sin(Omegak);Ygk = xk * Math.sin(Omegak) + yk * Math.cos(ik) * Math.cos(Omegak);Zgk = yk * Math.sin(ik);Rxf = -5.0 * Math.PI / 180.0D;Rzf = OmegaEdot * tk;Rx = new double[][]{{1D, 0D, 0D}, {0D, Math.cos(Rxf), Math.sin(Rxf)}, {0d, -Math.sin(Rxf), Math.cos(Rxf)}};Rz = new double[][]{{Math.cos(Rzf), Math.sin(Rzf), 0d}, {-Math.sin(Rzf), Math.cos(Rzf), 0d}, {0d, 0d, 1d}};Rgk = new double[][]{{Xgk}, {Ygk}, {Zgk}};XYZ = matrixMultiplication(matrixMultiplication(Rz, Rx), Rgk);Xk = XYZ[0][0];Yk = XYZ[0][1];Zk = XYZ[0][2];}//将计算结果添加到字符串outputDta.append(calculateHashMap.get("PRN") + "\t" + Xk + "\t" + Yk + "\t" + Zk + "\n");}//将字符串写入输出对象文件fileWriter.write(outputDta.toString());//回收资源bufferedReader.close();fileReader.close();fileWriter.close();}//计算儒略日static double calculateJD(double y, double mon, double d, double h, double miu, double s) {double day = d + (h * 3600 + miu * 60 + s) / 86400;double year = y + 4800;double m = mon;if (m < 2) {m = m + 12;year = year - 1;}double a = Math.floor(30.6 * (m + 1));double b = 0;if (y < 1582 || (y == 1582 && mon < 10) || (y == 1582 && mon == 10 && d < 5)) {b = -38;} else {b = Math.floor((y / 400) - (y / 100));}double c = Math.floor(365.25 * year);return (a + b + c + day);}//矩阵乘法static double[][] matrixMultiplication(double[][] A, double[][] B) {int line = A.length;        //矩阵A第一维长度(行)int column = B[0].length;   //矩阵B第二维长度(列)double[][] result = new double[line][column];   //定位结果矩阵大小for (int i = 0; i < A.length; i++) {            //矩阵A的行数for (int j = 0; j < B[0].length; j++) {     //矩阵B的列数double sum = 0;//矩阵的每一行乘以每一列for (int k = 0; k < B.length; k++) {sum += A[i][k] * B[k][j];}result[i][j] = sum;}}return result;}
}

五、实验结果

需要注意的是：

确保程序中使用的坐标系和单位的一致性。确保在计算过程中使用了正确的时间表达方式和相关的算法，以避免时间误差引起的计算错误。确保代码的可读性，采用良好的编码规范和注释。

基于卫星星历计算卫星在CGCS2000大地坐标系中的坐标相关推荐

基于MATLAB计算卫星位置
matlab卫星定位认识星历文件观测文件及位置计算星历文件观测文件 matlab面对对象卫星位置计算卫星高度角方位角计算(待更新...) 认识星历文件观测文件及位置计算初学者,欢迎指正批评. ...
RTKLIB学习总结（六）导航电文、卫星位置计算
文章目录一.导航电文 1.GNSS卫星信号的组成 2.导航电文的编排 3.遥测字(TLW) 4.交接字(HOW) 5.第一数据块 6.第二数据块 7.第三数据块二.卫星钟差钟漂改正 1.时钟校正参 ...
【python】读取卫星星历（RENIX 3.04）进行卫星位置的计算（北斗卫星专题）
最近的卫星导航数据处理,老师让我们进行卫星位置的计算,从而使用绘图工具进行对卫星星下点的轨迹进行绘图,这里首先的步骤是读取卫星星历数据,计算卫星位置. 这次的课程目标主要是针对北斗卫星,进行对卫星位置 ...
【Java】读取单颗北斗卫星导航星历文件进行卫星位置计算（RENIX3.04）
该文章是[python]读取卫星星历(RENIX 3.04)进行卫星位置的计算(北斗卫星专题)的一个其他语言实现的一个副本,其中可以计算出北斗GEO卫星的位置,进行检校. 其中文件的读取 ...
卫星过顶计算matlab,基于SGP4模型的卫星轨道计算.docx
基于 SGP4 模型的卫星轨道计算刁宁辉 ① , 基于 SGP4 模型的卫星轨道计算刁宁辉 ① ,刘建强 ① ,孙从容 ① ,孟鹏 ② (① 国家卫星海洋应用中心 ,北京 ...
实践：基于深度学习的卫星图像超分辨率处理技术（一）
本文由北邮@爱可可-爱生活老师推荐,阿里云云栖社区组织翻译. 以下为译文: Yann Lecun曾将无监督学习比作蛋糕,将监督学习比作蛋糕上的糖霜,声称我们仅懂得如何做糖霜却不知道怎样才能做出蛋糕. ...
计算卫星高度角、方位角
最小二乘定权.电离层对流层改正,都需要卫星的高度角.方位角.本章将介绍求解完卫星的地固坐标系的位置后,如何求解卫星的高度角.方位角. 卫星位置求解请参考之前的博客:卫星位置解算原理与程序设计参考书籍 ...
卫星位置计算小程序（C#版）
一.最终画面二.程序代码 using System; using System.Collections.Generic; using System.ComponentModel; using Sys ...
用matlab计算卫星的瞬时位置,用MATLAB计算GPS卫星位置.PDF
用MATLAB计算GPS卫星位置高新技术用MATLAB计算GPS卫星位置罗利娟杨乐 (西安翻译学院陕西西安 710061) 摘要:本文主要介绍了GPs测量数据的常用格式RINEX标准 ...

基于卫星星历计算卫星在CGCS2000大地坐标系中的坐标