一、北斗系统概述

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
地球（包含大气层）引力常数：	$\mu = 3.986004418\times 10^{14} m^{3}/s^{2}$
扁率：	f = 1/298.257222101
地球自转角速度：	$\Omega e = 7.2921150\times 10^{-5} rad/s$

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));

            //计算vk
            double 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);
            }

            //计算纬度幅角参数faik
            double faik = vk + Double.parseDouble(calculateHashMap.get("w"));

            //计算纬度幅角改正项duk、径向改正项drk、轨道倾角改正项dik
            double 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);

            //计算改正后的纬度幅角uk
            double uk = faik + duk;

            //计算改正后的径向rk
            double rk = A * (1 - e * Math.cos(Ek)) + drk;

            //计算改正后的轨道倾角ik
            double 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;
    }
}