卫星轨道计算

 

一． GPS观测量

接收机在观测相位和伪距数据的同时，还将广播星历和预报星历记录下来。接收GPS信号还能获取纳秒级精度的时间基准信号。

    由于接收机的型号很多，厂商设计的数据格式各不相同，国际上为了能统一使用不同接收机的数据，设计了一种与接收机无关的RINEX(The Receiver Independent Exchange Format)格式，目前已使用2号版本。下面分别介绍RINEX 2格式的广播星历文件、观测数据文件、和地面气象数据文件。

    RINEX 2格式的GPS数据文件的命名规则为：

     

其中： ～以4个字节表示的台站名；

      ～文件中第一组数据观测时间的年积日(例如：1月1日为001，2月2日为032)；

      ～该站该日收到的某类文件的顺序号，0表示只有一个；

      ～以两位数表示的年(例如：96表示1996年)；

      ～文件种类：

        ～观测数据文件；

        ～广播星历文件；

        ～地面气象数据文件。

    为了便于交流，RINEX 2格式的GPS数据文件均以①无带标；②ASCII码；③每个记录长度为80个字符，块大小为8000；录制在磁带上，磁带上的第一个文件是全部文件的目录。但目前国际上的IGS等组织是通过通讯方式(Internet网)，来快速地提取全球GPS长年观测站数据的，并将数据存在大型计算机中，使用着可通过Internet网任意提取。

    应注意，在RINEX 2格式的GPS数据中，时间均以GPST计，即与UTC要差一个整数跳秒数。

 

⒈ 广播星历文件

接收机锁定卫星并解出C/A码后，就能取得广播星历，即卫星坐标计算参数，在实时GPS应用中，它是必不可少的，大部分的工程网观测数据的后处理也采用广播星历。RINEX

 

2   格式的广播星历文件如下表2.1.1所示，作为例子，表中给出了PRN9和PRN17两颗卫星的广播星历数据，PRN表示GPS卫星的伪随机编号号码，GPS卫星在有些场合采用美国航空与航天局NASA(National Aeronautics and Space Administration)的编号。

表2.1.1 RINEX 2格式的广播星历文件

2              NAVIGATION DATA                         RINEX VERSION / TYPE EPHTORNX            Version 1.09        29-NOV-95 21:05     PGM / RUN BY / DATE      .1676D-07   .2235D-07  -.1192D-06  -.1192D-06          ION ALPHA      .1208D+06   .1310D-07  -.1310D+06  -.1966D+06          ION BETA      .133179128170D-06  .107469588780D-12   552960       39 DELTA-UTC: A0,A1,T,W     10                                                      LEAP SECONDS                                                             END OF HEADER  9 94 10 21  8  0  0.0-0.103851780295D-04-0.909494701773D-12 0.000000000000D+00     0.720000000000D+02 0.106062500000D+03 0.476841277575D-08 0.132076112444D+01     0.548548996449D-05 0.312971079256D-02 0.747293233871D-05 0.515371790504D+04     0.460800000000D+06 0.558793544769D-07-0.229012694900D+01-0.130385160446D-07     0.950477774712D+00 0.229593750000D+03-0.491558992251D+00-0.819034084998D-08     0.233938313166D-09 0.100000000000D+01 0.771000000000D+03 0.000000000000D+00     0.700000000000D+01 0.000000000000D+00 0.139698386192D-08 0.328000000000D+03     0.000000000000D+00 17 94 10 21  8  0  0.0-0.635907053947D-04-0.909494701773D-12 0.000000000000D+00     0.228000000000D+03 0.167187500000D+02 0.424946255961D-08 0.104717256943D+01     0.566244125366D-06 0.790101150051D-02 0.883266329765D-05 0.515369790649D+04     0.460800000000D+06 0.279396772385D-07 0.925235566518D+00-0.203028321266D-06     0.970438658460D+00 0.214781250000D+03 0.199825790573D+01-0.789747200969D-08     0.404659722397D-09 0.100000000000D+01 0.771000000000D+03 0.000000000000D+00     0.700000000000D+01 0.000000000000D+00 0.139698386192D-08 0.484000000000D+03     0.000000000000D+00

    上表中的前几行为表头，表头中的第60～80个字符是相应行的说明，如第一行的“2”是RINEX版本号，“NAVIGATION DATA”是指本文件类型为广播星历；第二行是生成该文件的单位、执行人、及文件形成日期。表头以“END OF HEADER”表示结束，表头结束有时也用空行表示。表头中注解为“ION ALPHA”和“ION BETA”的两行指相应的参数是电离层改正参数(采用差分模型相对定位时，该参数无用)，注解为“LEAP SECONDS”的行给出了GPST与UTC之间的跳秒数，在表头结束前可插入无限多的注解行，注解行的说明为“COMMENT”。

    表头结束后，每8行表示一颗卫星的广播星历，第一行的第1个数字是卫星的PRN号码(伪随机编号)。如表2.1.1中的第八行至第十五行是卫星PRN9的广播星历，对应量的含义如表2.1.2。

 

表2.1.2 RINEX 2格式广播星历数据含义

卫星PRN号，年，月，日，时，分，秒， ， ，

， ， ，

， ， ，

， ， ，

， ， ，

， ， ，

， ， ，

 

表中：

    ～卫星钟差常数项，时秒；

    ～卫星钟差漂移项，时秒/时秒；

    ～卫星钟差漂移速率项，时秒/时秒 ；

    ～平近点角的长期变化(近地点参数)，弧度/时秒；

    ～参考时刻的平近点角，弧度；

    ～扁心率；

    ～长半轴的平方根，米 ；

    ～参考时刻升交点赤经，弧度；

    ～参考时刻轨道倾角，弧度；

    ～近地点角距，弧度；

    ～升交点赤经在赤道平面中的长期变化(主要是由地球引力场的2阶带谐系数 ( )引起)，弧度/时秒；

    ～在星历参考时刻 在轨道延迹方向上周期改正余弦项的振幅，弧度；

    ～在星历参考时刻 在轨道延迹方向上周期改正正弦项的振幅，弧度；

    ～在星历参考时刻 在轨道径向方向上周期改正余弦项的振幅，米；

    ～在星历参考时刻 在轨道径向方向上周期改正正弦项的振幅，米；

    ～在星历参考时刻 轨道倾角(近似于法向)周期改正余弦项的振幅，弧度；

    ～在星历参考时刻 轨道倾角(近似于法向)周期改正正弦项的振幅，弧度；

    ～星历参考时刻(星期中的秒数)，秒；

    ～星历数据的年龄；

    ～轨道倾角变化率，弧度/时秒；

    ～ 上存在哪些码指示；

    ～GPS星期数；

    ～ 上P码伪距指示；

    ～本广播星历精度指标，米；

    ～卫星是否健康指标；

    ～电离层群延迟改正参数，时秒；

    ～卫星钟数据年龄;

    ～信息传送时间(与接收机对接收到的卫星信号解码有关)，时秒。

 

 

⒉ 观测文件

    接收机测得的相位和伪距观测值均记录在观测文件中，下表2.1.3是例子文件，共有2个历元的观测数据。

 

表2.1.3 RINEX 2格式的观测数据文件

2              OBSERVATION DATA    GPS                 RINEX VERSION / TYPE TB2RNX              xxxxxxxxxxxx         95-07-20 22:20:20  PGM / RUN BY / DATE Turbo SII rinex formatter Version: 95.5.19                  COMMENT MODE : STATIC                                               COMMENT                                                             COMMENT 1                                                           MARKER NAME WAN                                                         OBSERVER / AGENCY 153575902           TURBO SII           Production unit     REC # / TYPE / VERS 153575902           TURBO SII                               ANT # / TYPE  -2852389.2261  4650364.5453  3293350.3434                  APPROX POSITION XYZ         0.1120        0.0000        0.0000                  ANTENNA: DELTA H/E/N Original slant height(m) :    0.1120                        COMMENT      1     1                                                WAVELENGTH FACT L1/2      5    C1    L1    L2    D1    P2                        # / TYPES OF OBSERV                                                             END OF HEADER  95  7 19  0 21  0.0000000  0  5 27 28 15 31 19   21805891.33516        -0.28716         0.07713         0.00000  21805895.08513   23733969.59714        -0.19314         0.66912         0.00000  23733973.16712   20746137.25719         0.01719         0.54216         0.00000  20746140.99416   20903455.21318        -0.28018         0.99414         0.00000  20903458.41214   20681238.04917        -0.32817        -0.25610         0.00000  20681238.04910  95  7 19  0 21 30.0000000  0  5 27 28 15 31 19   21794932.63606    -57587.90306    -44873.38602         0.00000  21794935.76502   23750223.73904     85421.79604     66563.21801         0.00000  23750225.60401   20736273.69909    -51834.40309    -40389.90405         0.00000  20736277.71705   20913805.23508     54389.67808     42382.77504         0.00000  20913808.93304   20682314.76407      4906.19307      3823.00700         0.00000  20682314.76400

 

    与广播星历表一样，观测文件也有一个表头，表头也以“END OF HEADER”或空行表示结束，表头中每行的60～80字符为给出本行内容的说明，如表2.1.3所示，观测文件的表头中列出了RINEX版本号、形成文件的单位、人员、点号、点名、观测者与单位、接收机号码类型版本、天线号码类型、台站近似坐标、天线L1相位中心与点位的关系(向上、向东、向北的偏离量)、采样间隔、波长因子、观测值种类数及观测值类型、第一个观测历元时刻、最后一个观测历元时刻、测得卫星的号码，说明为“COMMENT”的注解行等。

    其中的波长因子为1指采用相关技术恢复载波，整周模糊度和失周数只能是整数；波长因子为2指采用平方技术生成载波，整周模糊度和失周数可能是0.5周的倍数。观测类型中的：

       L1～L1上的载波相位；

       L2～L2上的载波相位；

       C1～L1上的C/A码伪距；

       P1～L1上的P码伪距；

       P2～L2上的P码伪距；

       D1～L1上的频率变化；

       D2～L2上的频率变化。

 

    表头结束后，是观测各历元的观测数据，每历元数据由一个历元时间行和多个观测数据行组成，每个历元时间行的内容为：

年，月，日，时，分，秒，质量标记，卫星数，卫星号码，卫星号码，... ，钟差

其中：

    ①卫星数指该历元观测到的卫星总数，紧接着的是观测到的卫星号码序列。

    当质量标记(Event Flag)为：

      ～该历元观测数据正常；

      ～在前一历元和本历元之间停电；

      2～开始移动天线；

      ～在新的台站重新开始观测(紧接着会出现新的点号)；

      ～以下会出现新的头信息；

      ～其它外部事件；

      ～出现失周指示。

    ②卫星号码计为 ：

      ：卫星系统：

          ～GPS；

          ～GLONASS系统(前苏联研制的类似于GPS的定位系统)；

          ～多普勒卫星定位系统。  

      ：如果是GPS系统，PRN编号；

           如果是GLONASS系统，通道号；

           如果其它卫星系统，两位数的编号。

    ③历元行中的钟差(选项)，是接收机钟差，处在68～80位，如果此项存在，则应对历元时刻、测得的伪距相位作如下修正：

      历元时刻 = 给出历元时刻 - 钟差

      伪距 = 测得伪距 - 钟差×光速

      相位 = 测得相位 - 钟差×频率

 

    历元行后的几行(行数等于卫星总数)，是对应的每个卫星观测值，每行中的数据个数等于表头中的观测值种类数，数据类型按表头中的观测值类型排列，每个观测值后面紧跟着一个一位数的失周指示和一个一位数的信号强度指示。相位观测值的单位是周，伪距的单位是米，当某个数据没有测到时，计为0.0或空格。

    当失周指示(取值范围是0～7)为 或空格时数据正常。

    当信号强度指示为：

      或空格～正常；

      ～信号最弱情况；

      ～最理想的信噪比；

      ～信号最强情况。

 

⒊ 气象文件

    在高精度观测时，有时需要记录气象数据，RINEX 2格式的气象数据文件如下表2.1.4所示：

 

表2.1.4 RINEX 2格式的气象数据文件

2              METEOROLOGICAL DATA                     RINEX VERSION / TYPE TB2RNX              xxxxxxxxxxxx         95-07-20 22:20:20  PGM / RUN BY / DATE 1                                                           MARKER NAME      3    PR    TD    HR                                    # / TYPES OF OBSERV                                                             END OF HEADER  95  7 19  0 10 00  987.1   10.6   89.5  95  7 19  1 10 00  987.2   10.9   90.0  95  7 19  2 10 00  987.1   11.5   89.0

 

    表头及结束符与星历和观测数据文件类似，记录的是各气象观测时刻(年月日时分秒)，台站周围的气压(毫巴)、干温(摄氏度)、相对湿度(百分比)。

 

二． 由广播星历计算卫星位置

    如果要计算时刻 某卫星的空间坐标，根据前面的广播星历格式，读出该卫星广播星历的有关参数，按如下步骤计算：

1. 求长半轴

      

2. 计算平角速度

      

式中的 为地球引力常数。

3. 计算从需要时刻到参考时刻的时间差

      

4. 改正平角速度

      

5. 计算平近点角

      

6. 按下式迭代计算偏近点角

      

7. 由下两式计算真近点角

      

8. 计算纬度参数

      

9. 周期改正项

      

10. 计算改正后的纬度参数  

      

11. 计算改正后的向经

      

12. 计算改正后的倾角

      

13. 计算卫星在轨道平面内的坐标

      

14. 改正升交点的经度

      

式中的 是地球自转的角速度。

15. 最后计算卫星在地固系中的坐标

      

    根据以上过程计算卫星位置时，应注意，求出卫星坐标表示在地固系中，计算时刻 必须是以GPST计的卫星信号发射时刻，它由信号接收时刻减去传播时间延迟(需迭代求得，见后)后得到。当 大于302400秒时，因在 中减去一个GPS星期相应的秒数604800秒，当 小于0秒时，应在 中加上604800秒。

通常认为按以上过程求出的坐标是表示在WGS84中的，但由于其误差很大，个别情况甚至超过100米，因此可认为该坐标表示在任一地固参考系CTS中。

 

16. 计算GPS信号从卫星发射时刻的时间（减去到接收机传播时间的延迟）

    GPS观测值是卫星和接收机振荡器频率及传播时间延迟 的函数，传播时间延迟包括几何传播时间延迟 和大气传播时间延迟 ，即：

                                                          2.2.15

    我们进行GPS测量需要求定的参数均包含在 中，下面讲述 的求定。

    严格来说， 的求定应该在广义相对论框架下的以太阳系质心为原点的坐标系中计算。但实用时，在地心惯性坐标系下计算就能满足精度要求，等价地，也可在地固坐标系中计算。

    在定位应用中，卫星轨道是已知的，台站坐标是未知的，但近似值是知道的(即使没有其它已知信息，也能从观测文件中取得台站的近似坐标)，数据处理的迭代过程就是将台站坐标从近似值算到准确值。在定位定轨同时进行时，卫星轨道的近似值也是知道的，在没有其它已知信息的情况下，可从广播星历文件中取得。

    如果在地固坐标系中计算在 时刻台站 收到的卫星 的信号的几何传播时间延迟为 ，按电磁波传播定律：

                                             

式中的 是卫星在信号发射时刻 时，表示在地固坐标系中的坐标； 是台站在地固坐标系中的坐标； 是电磁波在真空中的传播速度(光速，299792458米/秒)。由于包含在2.2.16式中的信号发射时刻 是未知的，要在计算出 之后才能求出，因此需要一个迭代过程才能求得 ，具体步骤是：

    ①按第 次迭代求得的几何传播时间延迟值 (在第一次迭代时，可假定  秒，因为卫星离地面的高度近似为20000公里)，求出卫星信号发射时刻 ：

        

      

式中的大气传播时间延迟 可按大气改正模型计算(详见后)。

    ②按卫星星历表内插求得 时刻卫星坐标 。

    ③按上式计算 ：

     

    ④当 与 之差大于限差时，回到①，当 与 之差小于限差时，则认为已经收敛，并以 的值代替 ，加上 即得 ，该限差应取为小于等于 的值，因为它对相位的影响是它与频率的乘积，而频率的量级是 ，GPS接收对相位的量测精度可达0.01周。

以上的传播时间延迟是在地固系中计算的，也可以在惯性系中计算，可按方便选用，如果星历是表示在地固系中的(如：广播星历和NGS格式的精密星历)，而又不进行定轨，则选择在地固系中计算较为方便；如果要求定轨道或星历已经表示在惯性系中，则在惯性系中计算较为方便，计算时需要在每个历元，将台站坐标旋转到惯性系中，而避免计算量很大的将星历从惯性系旋转到地固系的计算。

 

17. 中心大气改正

    从地表至离地面80公里高的这一层大气中，原子和分子处于中性状态，故称为中性大气，也可称为对流层，有的文献将它又分为对流层和平流层两层，但由于GPS观测仅采用高度角大于15度的观测值，分一层和分两层在效果上并无区别。中性大气层使得电磁波的传播时间增加，称为中性大气延迟，在天顶方向就可达到2.5米左右，并随高度角的增加而加大。中性大气延迟分两个部分，由大气中所有大气分子的偏振位移引起的称为干项，由水分子的偶极距引起的称为湿项。其中干项比较稳定，用合适的模型可得到较好的改正，湿项引起的附加延迟要小得多，只有几十厘米，但其变化很不规则，没有高精度的改正模型，采可昂贵的水气辐射机可以测定，大部分情况下采用引入待估参数的方法来处理。

    天顶方向大气延迟改正 可采用SAASTAMOINEN 1973、MARINI、CHAO等模型来计算，一般方向上的改正还必须乘以一个以天顶距 为变量的映射函数 ，映射函数有Lanyi、CFA等模型。在GPS数据处理中，观测数据量非常大，我们一般只采用高度角大于15 的观测数据，所以大部分的现有天顶改正和映射函数模型均可采用。经常采用的SAASTAMOINEN 1973模型如下：

                                   

其中干项天顶延迟：

       

湿项天顶延迟：

       

     

     

    由于4.1式中的湿项模型不准，一般可引进一个待定参数 ，使4.1式变为：

                             4.1.2     

                                   

上式中的：

     

       

     

       

     

以上的 为地面气压(毫巴)， 为地面温度(摄氏)， 为水气压(毫巴)， 为相对湿度， 为台站的地心纬度， 为台站的大地水准面高(公里)。

    水汽辐射计可以用以精确测定大气层中积累的水蒸汽和云雾对电磁波传播路径增长的影响，模型如下：

    以 表示水汽辐射计测得的任意方向上的天空亮温度，则大气湿项引起的电磁波传播路径增长 (即前面的 与光速的乘积)可由以下一组公式计算：

                 

上式中的 和 为两个波段的频率， 和 为对应的亮温度， 可近似地取为 ， 是宇宙背景的辐射温度， 是传播路径对应的权函数， 和 为水蒸汽和氧的吸收系数， 和 可用一个标准的大气模型以及地面所测温度来确定，近似地，  ( 为地面温度， )。水汽辐射计测得的 的精度可达到1.5毫米左右。

 

18. 电离层改正

    地球表面60～1000公里的这层大气，由于太阳辐射，其中的原子被电离成大量的正离子和电子，构成电离层。电磁波通过电离层产生的时间延迟 为：

                                                     4.1.5     

式中的 为观测频率， 为光速， 为传播路径上的电子总含量， 为经典电子半径。

    由于4.1.5式中的 与太阳辐射压有关，昼夜可相差一个数量级，很难用模型来精确表示。但可以看出它产生的传播时间延迟与频率的平方成反比，故可以用双频观测来消除其影响。

    如假设 、 、 、 为某GPS接收机在同一历元测得的 和 两个波段上的伪距和相位，则电离层对它们的影响分别为 、 、 、 ，这里的 对一个接收机对同一卫星在同一历元的几种观测量是常数， 和 为两个波段的频率。因此在长距离GPS相对定位中，通常将两个波段的相位观测值线性组合后的观测量 作为观测值。

                                                    

式中的 ，即L1与L2的频率之比。

 

三．伪距差分定位模型

    基准站R测得至GPS卫星J的伪距为

   

其中 为基准站R至第J颗卫星之间的真实距离； 是GPS卫星星历误差所引起的距离偏差； 为接收机时钟相对于GPS时间系统的偏差； 是J颗卫星时钟相对于GPS时间系统的偏差； 电离层时延所引起的距离偏差； 是对流层时延所引起的距离偏差；c是为电磁波的传播速度。

根据基准站的已知坐标和GPS卫星星历，可以精确算得真实距离 ，而伪距 是基准站接收机测得的，则伪距的改正值

在基准接收机进行伪距测量的同时，流动站接收机K也对第J颗卫星进行了伪距测量，流动站接收机所测得的伪距

我们将基准站所测得的伪距改正值带入上式（即上两式相加）得：

当流动站与基准站相距在一定距离范围之内，我们可以近似认为：

 

则上式变为：

(6)式中共有四个未知数，它们分别是流动站K的三维坐标 和GPS接收机钟差引起的改正项 。如果基准站和流动站某历元共视四颗以上卫星，就可以根据上式建立误差方程式：

  

式中（ ）是流动站K的三维近似坐标；（ ）是卫星J发射信号时刻的三维坐标，其可根据卫星星历计算而得到； ；j=1,2,….n  n= 。

然后按最小二乘法求解流动站该历元的三维坐标 和三维坐标中误差 。