GNSS地面基准站数量、测站分布和观测质量是影响卫星定轨精度的主要因素。文献[1]全面地分析了基准站分布对GLONASS卫星定轨与钟差精度的影响。研究表明，基准站布局对区域导航系统具有重要作用，证实了测站的分布和几何强度对定轨精度的影响^[2-4]。目前，IGS(international GNSS service)全球跟踪站达到500多个，如何利用现有站点优化选取来提高定轨精度，控制计算成本，进而提高GNSS实时产品服务能力，还有待更深入研究。

当前，普遍采用格网法进行基准站全球均匀选站^[5]。其优点是可相对简单直观地获取测站的均匀分布构型，以达到较高定轨精度，然而，当同时考虑多个影响因素(如几何分布、站点稳定性等)时，很难实现全局最优化设计，不能保证几何构型达到最佳^[6-7]。此外，该方法属于人机交互操作模式，人为因素干涉大，比较费时。且传统格网法缺少定轨地面构型的评价指标。有研究提出用几何精度衰减因子(geometric dilution of precision, GDOP)值来反映测站分布均匀程度^[8-10]。GDOP值越小，测站分布越均匀，几何构型越好。文献[5]研究讨论了卫星定轨地面站几何构型对定轨精度的影响，地面站点的MDOP (定轨动力学参数精度衰减因子)越小，卫星定轨精度越高。MDOP值与地面站的均匀分布程度密切相关，且直接影响卫星轨道精度。

本文提出基于格网控制概率分配的测站随机优化选取方法，实现定轨测站构型自动优化筛选。在兼顾测站多个质量因子和测站几何分布的基础上，引入基于地面站均匀分布的评价指标来衡量测站构型的优劣。该方法可在有限时间内获取更好的地面测站分布构型，提高定轨参数、钟差参数和地球自转参数等的解算精度；能够快速自动选取高质量和最佳几何分布的测站。

1 测站分布对定轨的影响

当利用地面观测站进行GNSS定轨时，卫星定轨的几何观测距方程可简写为：

$ \rho_{j_{k}}^{i}=R_{j_{k}}^{i}+f(p)+\varepsilon_{j_{k}}^{i} $

(1)

式中，ρ_jkⁱ为第k个历元卫星i与测站j之间实际距离；$R_{j_{k}}^{i}=\left[\left(x_{k}^{i}-X_{j k}\right)^{2}+\left(y_{k}^{i}-Y_{j_{k}}\right)^{2}+\left(z_{k}^{i}-\right.\right.Z_{j k} )^{2} ]^{\frac{1}{2}}$表示测站与卫星之间的几何距离；(x, y, z)表示卫星坐标；(X, Y, Z)表示测站坐标; ε_k表示观测噪声；f(p)表示钟差、模糊度、对流层、地球自转参数等附加模型参数。下文仅讨论轨道位置参数有关的几何图形矩阵，即假设这些参数均已解算得到。将式(1)按照一阶泰勒级数展开，得到误差方程式的线性化形式为：

$ \Delta \boldsymbol{\rho}_{k}=\boldsymbol{H}_{k} \Delta \boldsymbol{X}+\boldsymbol{\varepsilon}_{k}, k=1,2 \cdots N $

(2)

式中，H_k为第k个观测历元对应的几何图形矩阵，表达式为：

式中，n表示测站数；k表示观测历元。对式(2)按照最小二乘准则对卫星位置的修正值求解得^[11]：

$ \Delta \boldsymbol{X}=\left(\boldsymbol{H}^{\mathrm{T}} \boldsymbol{H}\right)^{-1} \boldsymbol{H}^{\mathrm{T}} \Delta \boldsymbol{\rho}_{k} $

(3)

其中，$\boldsymbol{H}=\sum\limits_{k=1}^{N} \boldsymbol{H}_{k}$, 由式(3)的误差引起的定轨位置误差为^[12]：

$ \sigma_{r}=\sqrt{\sigma_{x}^{2}+\sigma_{y}^{2}+\sigma_{z}^{2}}=\sigma_{0} \sqrt{\operatorname{tr}\left(\boldsymbol{H}^{\mathrm{T}} \boldsymbol{H}\right)^{-1}} $

(4)

$ \mathrm{MDOP}=\left[\operatorname{tr}\left(\boldsymbol{H}^{\mathrm{T}} \boldsymbol{H}\right)^{-1}\right]^{\frac{1}{2}} $

(5)

式中，tr表示矩阵对角线上的元素和，即矩阵的迹；σ₀为距离测量中误差。当中误差一定时，MDOP的大小主要取决于H矩阵结构，而矩阵H由卫星与基准站之间的方向余弦构成，因此矩阵H的结构与测站相对于卫星的空间分布有着密切的联系。

地面站对卫星定轨的几何构型随着时间不断变化(见图 1)。在卫星定轨中，常使用一定弧段的观测资料积累实现定轨，此时，多个地面站、多个历元对卫星的累积观测资料给出了卫星定轨的几何信息。对应第j个测站，其几何图形矩阵可参数化为：

$ \begin{array}{c}{\boldsymbol{H}_{k, j}=\left[\sin \theta_{j} \cos \varphi_{i} \quad \sin \theta_{j} \sin \varphi_{i} \quad \cos \theta_{j}\right]} \\ {j=1,2 \cdots n}\end{array} $

(6)

式中，θ_j为卫星与地心连线和卫星与测站分布组成的圆锥半角；φ_i为卫星与测站分布的连线投影到地面与地球直角坐标系x轴的夹角。研究表明^[13]，假设当地面测站均匀分布(如构成嵌套圆锥构型)时，可得：

$ \begin{array} [c]{c} \left(\boldsymbol{H}_{k}^{\mathrm{T}} \boldsymbol{H}_{k}\right)^{-1}=\\\operatorname{diag}\left(\frac{1}{\sum\limits_{c=1}^{T} \frac{n_{c}}{2} \sin ^{2} \theta_{j}}, \frac{1}{\sum\limits_{c=1}^{T} \frac{n_{c}}{2} \sin ^{2} \theta_{j}}, \frac{1}{\sum\limits_{c=1}^{T} n_{c} \cos ^{2} \theta_{j}}\right) \end{array} $

(7)

式中，n_c为第c个圆锥上的测站数；T为圆锥面总数。

当地面测站均匀分布在地球表面时，随着卫星的运动，其可见地面站都均匀可见，在此理想情形下，可累加得到：

$ \left(\boldsymbol{H}^{\mathrm{T}} \boldsymbol{H}\right)^{-1}=\left[\sum\limits_{k=1}^{\mathrm{N}}\left(\boldsymbol{H}_{k}^{\mathrm{T}} \boldsymbol{H}_{k}\right)\right]^{-1} $

(8)

由式(5)和式(8)可知，观测弧段越长，即观测历元越多，MDOP的值越小。将式(5)进一步展开得：

$ \begin{aligned} \mathrm{MDOP}^{2} &=2\left(\sum\limits_{k=1}^{N} \sum\limits_{c=1}^{T} \frac{n_{c}}{2} \sin ^{2} \theta_{j}\right)^{-1}+\\ &\left(\sum\limits_{k=1}^{N} \sum\limits_{c=1}^{T} n_{c} \cos ^{2} \theta_{j}\right)^{-1} \end{aligned} $

(9)

根据式(9)可知，在一定观测弧段，当测站数和圆锥面总数相同时，θ_j变大，即测站的分布较开阔时，MDOP值会变小，定轨精度提高。当测站分布均匀时，圆锥面总数会增多，T变大，MDOP的函数值减小。观测历元N越大，测站数n_c越多，都会使得MDOP的值变小。

可见，随着观测时间的增长、测站数目的增多，定轨的几何图形在不断改善，即定轨精度不断提高。卫星可见性和数量一定的条件下，测站分布越广，定轨效果越好。因此考虑测站分布和质量，完善测站的选取方法，对定轨精度的提高有重要影响。

2 全球测站均匀选取方法 2.1 格网控制法

目前常采用选站的方法是格网法，其基本思想是采用经纬度格网划分。该方法主要是根据选取的测站数量，用合适的经纬度格网划分全球，把多个离散站点划分在不同的区域。经纬度的格网数num用式(10)计算：

$ \operatorname{num}=\left[\frac{360}{10 n} \times \frac{180}{10 n}\right], \quad n=1,2 \cdots n $

(10)

式中, [ ]表示向上取整。

由于测站在全球分布极其不均匀，格网中站点数目也不同，这时一般需要综合考虑每个测站位置和质量，选取每个格网中较为理想的定轨测站构型。需要指出，格网法实际上很难实现一些素数问题的划分。

该方法主要考虑测站分布，很难同时权衡考虑测站稳定性和测站观测质量等信息，人为控制因素较多，具有一定的人为性。针对这种弊端，本文基于组合优化设计思想，针对全组合的候选构型数目巨大的问题，提出了导航卫星定轨测站构型的随机组合优化方法。

2.2 格网控制概率下的测站随机优选方法

基于以上问题，本文提出格网控制概率下的测站筛选方法，把概率统计思想引入测站选取中，即在格网法基础上，兼顾测站数据质量、稳定性和地理分布等信息，为每个测站综合分配一定的概率。基本思想如下：

1) 测站概率的格网控制

全球基准站的分布有地区性差异，为避免选取的测站集中在某一区域，需对测站位置进行整体控制，提高全局解随机优化收敛能力。将全球按照式(11)的格网均匀分区，统计格网中测站的数量及质量等信息。

$ J W=180 / \sqrt{n / 2} $

(11)

式中，JW表示经纬度网格的大小，单位为(°)；n为测站数。

2) 测站概率二次分配

当站点精度高，观测质量好，测站分配较大概率，概率分配计算公式为：

$ G_{j}^{l}=\frac{b {\sigma_{0 j}}^{2}}{\sum\limits_{j=1}^{n_{l}} {\sigma_{0 j}}^{2}}+\frac{d m p_{j}}{\sum\limits_{j=1}^{n_{l}} m p_{j}} $

(12)

式中，σ_0j²=σ_X²+σ_Y²+σ_Z², 其中σ_X、σ_Y、σ_Z为站点坐标分别在X、Y、Z方向的中误差；G_j^l表示第l个格网中j点的概率；n_l为格网中的点数; mp_j为j点的多路径误差值；b、d表示经验值。

3) 蒙特卡洛随机实验

按照随机选取的原则，首先给每一个测站分配相应的概率：

$ F_{j}=G_{j}^{l} \frac{1}{\mathrm{num}_{-} \mathrm{block}} $

(13)

式中，num_block为含有测站的格网数；F_j为测站j的概率。

将所有测站作为实验总体，样本为S，从中随机选取n个测站，即随机选取一组测站列表：

$ {\rm select}_{-} \text { list }=\operatorname{randsrc}(n, 1,[S ; F]) $

(14)

式中，F为包含所有测站概率的矩阵；randsrc为按照概率F大小在总样本S选取n个测站的算法；select_list为测站列表。

4) 实验样本的筛选指标

当控制点均匀分布时，控制图形几何中心位置的GDOP可达到极小值^[14]，为此可以计算地面n个测站对地心的GDOP值作为衡量地表测站均匀分布程度的指标。GDOP值越小，测站分布越均匀，同时几何图形对称性越好，此时定轨构型越好，每颗卫星定轨的动力学几何因子MDOP值越小。需要指出，构型优劣的评价指标并不唯一，例如，使用构型的体积最大指标、测站网间平均距离最短等指标。

本文提出加权GDOP准则，加权GDOP越小，选择测站的组合越好^[7]，即得到测站列表。这种方法不仅能考虑到测站的质量和几何构型，而且还能使两者达到很好组合，使得定轨结果达到最佳。加权GDOP计算如下：

$ \boldsymbol{A}=\left[\begin{array}{cccc}{X_{1}} & {Y_{1}} & {Z_{1}} & {1} \\ {X_{2}} & {Y_{2}} & {Z_{2}} & {1} \\ {\vdots} & {\vdots} & {\vdots} & {\vdots} \\ {X_{n}} & {Y_{n}} & {Z_{n}} & {1}\end{array}\right] $

(15)

$ \begin{array} [c]{c} \boldsymbol{P}=\operatorname{diag}\left(\frac{b s_{1}}{{\sigma_{0,1}}^{2}}+\frac{d s_{2}}{m p_{1}}, \frac{b s_{1}}{{\sigma_{0,2}}^{2}}+\frac{d s_{2}}{m p_{2}}\right) \cdot\\ \frac{b s_{1}}{{\sigma_{0, n}}^{2}}+\frac{d s_{2}}{m p_{n}} \end{array} $

(16)

$ \mathrm{GDOP}=\operatorname{tr}\left(\boldsymbol{A}^{\mathrm{T}} \boldsymbol{P} \boldsymbol{A}\right)^{-1} $

(17)

式中，s₁和s₂分别表示所有测站方差和多路径误差的中位数; b、d为经验值; P表示每组测站的权重。依据抗差估计理论^[15]，当bs₁/σ_{0, n}²+ds₂/mp_n大于1时直接赋权重为1，小于1时按照实际值输出。避免有些点对几何构型改善明显，但质量较差，测站给的权重太小，从而使加权GDOP值过大，丧失该测站在定轨中应有的价值，进而降低该指标的有效性^[16]。

3 实验分析 3.1 算法实现与流程设计

本文选取的实验数据来源于IGS网站，共203个站，测站精确坐标从SINEX文件获得。基于格网控制概率的随机优化方法，设定实验的采样数(本文选择10 000次)，然后按测站加权GDOP值大小排序，选出最小加权GDOP值的测站组合，即要选取的测站。其流程如图 2所示。

3.2 数据处理与实验分析

选择处理时间段为2016年年积日为180~194共15 d GPS观测数据。采取以上方法选取测站后，以IGS在官方发布的精密卫星的最终轨道作为参考评价的指标^[17]，数据处理策略见表 1。

按照随机优选方法进行10 000次采样，每组测站实际最小加权GDOP值和理论值与测站数量的关系如图 3所示。

从图 3可以看出，总体上，实际GDOP值与理论GDOP值随着测站增加而逐渐接近，测站数目达到一定数量，两者的变化曲线最终重合。随着测站数增加，加权GDOP的值逐渐减小，测站构成的几何图形逐渐变好。当测站数比较少，特别是小于60个站时，GDOP值的变化率比较大，说明在一定范围内，测站数增加，几何构型改善比较明显。当测站数达到60个以上，加权GDOP值基本不再发生变化，它的变化率趋向于0，这时再增加测站对几何结构改善随着测站数目增多变化缓慢。

以30个站为例，最小GDOP测站分布与卫星运行的轨迹如图 4所示(其中蓝色三角为测站，红色轨迹为卫星G01~G10运行轨迹)。选取样本中10个加权GDOP最小和10个最大的测站列表进行实验，以便对结果做统计分析。选取10、20、…、90个测站对GPS卫星进行精密定轨，计算轨道结果与IGS最终轨道比较，限于篇幅，给出部分精度图。

图 5中，横坐标表示实验的组数，1~10和11~20分别表示选取最小和最大的加权GDOP值的实验组各10个，纵坐标表示GPS卫星定轨精度结果的均值。可以看出，当测站数相同时，前10组实验的定轨结果明显好于后10组，加权GDOP值的大小能够在一定程度上表示定轨精度的好坏。测站数一定，加权GDOP值越小，定轨精度越好。其中当测站数为10个时，最小和最大加权GDOP实验组的平均定轨精度分别为10.53 cm、43.14 cm，两者精度相差32.61 cm。测站数为30时，两者的平均定轨精度分别为2.13 cm、3.57 cm；90个测站时，两者的定轨精度分别为1.10 cm、1.32 cm，仅差0.22 cm。说明当测站达到一定数量，几何构型对导航卫星定轨精度的影响逐渐减弱。前10组计算的结果波动小，后10组相对波动较大。可以得出，加权GDOP值越小，计算结果越稳定。图 6表示加权GDOP值与定轨精度的关系图，可以看出定轨精度与加权GDOP的值成负相关。

随着测站数增多，GPS定轨精度会逐渐升高。选取30个测站时，GPS精密定轨的精度为2.13 cm。测站数小于60时，增加测站，定轨精度提高显著；测站数多于60时，定轨精度不会再有较大的提升。

分别以30和60个测站为例，计算15 d GPS精密定轨的结果，用传统格网法与本文提出方法轨道精度进行对比分析，如图 7所示。可以看出，随机优化选站的方法精度整体要高于格网法，在30个测站时，格网法平均定轨精度为3.64 cm，随机优化算法的精度为2.15 cm，提高了1.49 cm。60个测站时，格网法的精度为1.38 cm，本文方法的精度为1.26 cm，整体精度提高9.5%。可以看出，测站数量越少，本文方法的优势越明显。

用随机优化的方法选取90、160和203个测站计算的轨道精度如图 8所示。

从图 8中可以看出，测站数为90、160和203时，定轨精度与测站数为60时相当，定轨精度并没有随测站的增加而提高，测站数为160时，平均轨道精度为1.03 cm，全部203个测站参与解算时，精度为1.08 cm。可见，当测站数达到一定程度，测站数增加并不会显著提高定轨精度，反而可能会使定轨精度降低。而且随着测站数增多，需要解算的定轨参数会大大增加，解算时间会加倍增长。

4 结语

本文提出的基于格网控制概率下的随机优选方法，能够快速自动选取几何分布与测站质量占优的测站列表，可广泛应用于实时精密轨道和预报轨道的解算中。通过该方法统计分析了15 d观测数据的定轨结果，进一步得出以下结论：①加权GDOP值越小，测站与卫星组成的几何构型越好，定轨的精度越高。②在一定范围内，随着测站数的增加，选取的测站数大于60时，定轨几何结构信息改善变化缓慢。③选取60个测站，定轨的精度能达到1.26 cm，与基于所选取的203个测站的定轨精度相当。并且当选取160个测站时，其精度比全部测站参与解算的精度要高。表明测站数目多，解算的轨道精度不一定高。

该方法应用效果优于格网法，特别是选取少量测站时，例如30个测站，定轨精度提高显著。结合神经网络法或蚁群，有望进一步提高该方法全局最优解的搜索效率。