提高导航卫星精密定轨解算效率是GNSS卫星精密定轨技术的研究重点之一^[1-2]。Ge等^[3]采用消参方法，在每个历元时刻消去法方程中的过时参数，极大地减少了最终法方程解算参数的个数，从而有效提高了导航卫星精密定轨解算的效率，该方法是从算法层面对解算效率进行优化。随着CPU的迭代发展，部分学者开展基于CPU技术提高计算效率的研究，包括基于CPU的多核、多线程并行运算技术^[4-5]和基于CPU分层存储架构的方法^[6]。

目前，关于提高GNSS数据处理软件计算效率的研究几乎都是基于X86架构的Intel CPU，涉及国产ARM架构CPU^[7]上效率提升的实验和研究很少。在科学计算中，线性代数库是运行效率优化的重要工具，Lapack能以高效的方式实现科学计算的接口规范^[8]，并且适用于X86、ARM和PowerPC等不同架构平台。Intel CPU发展至今已经很成熟，Intel可实现适应其CPU的高效运算库Intel Math Kernel Library(MKL)。表 1为导航卫星精密定轨处理软件分别运行于Intel CPU中性能突出的X86架构Intel(R) Xeon(R) Gold 6134 CPU(表中Intel)和国产CPU中性能突出的飞腾ARM FT-2000+CPU(表中ARM)的运行耗时，软件运行时采用相同的数据和场景。可以看出，采用Lapack时，Intel CPU的解算效率是ARM CPU的3倍多；采用MKL后，Intel CPU的解算效率是ARM CPU的5倍多。由此可知，ARM CPU还需要寻求更合适的优化方法以提升导航卫星精密定轨的解算效率。

本文选取国产飞腾CPU的最新型号FT-2000+进行卫星精密定轨解算效率优化研究。

1 GNSS卫星精密定轨解算 1.1 解算原理

设历元i的第j条观测是接收机r与卫星s的伪距相位观测，添加动力学约束的观测方程为：

$\begin{gathered} P_j\left(t_i\right)=f\left(\vec{X}_{\mathrm{s}}, \vec{X}_{\mathrm{r}}, \overrightarrow{\mathrm{ERP}}\right)+\mathrm{ZTD}_{\mathrm{r}} \cdot \operatorname{map}+ \\ \mathrm{AMB}_{\mathrm{r}}^{\mathrm{s}}+\mathrm{clk}_{\mathrm{r}}\left(t_i\right)+\mathrm{clk}_{\mathrm{s}}\left(t_i\right)+\varepsilon \end{gathered}$

(1)

式中，P_j(t_i)为相位值的双频无电离层组合(伪距与之相同，不同的是右方程无模糊度参数)；$\vec{X}_{\mathrm{s}}$为卫星s的待估动力学参数，包括卫星初轨参数、太阳辐射压参数，其中太阳辐射压参数为ECOM1模型5参数；$\vec{X}_{\mathrm{r}}$为接收机r的位置参数；$\overrightarrow{\mathrm{ERP}}$为待估的地球定向参数(包括极移X_p、Y_p，极移速率$\dot{X}_p、\dot{Y}_p$和日长变化LOD等)；ZTD_r为接收机r的天顶方向对流层延迟；map为映射函数；AMB_r^s为站星对的模糊度参数；clk_r(t_i)和clk_s(t_i)分别为接收机r和卫星s的钟差参数；ε为测量随机误差。由于卫星钟差和接收机钟差具有随机性，因此将二者作为白噪声进行逐历元估计，clk_r(t_i)、clk_s(t_i)的有效期仅为当前历元。

将式(1)在近似值处线性化得到：

$\begin{aligned} P_j\left(t_i\right)= & f\left(\Delta \vec{X}_{\mathrm{s}}, \Delta \vec{X}_{\mathrm{r}}, \Delta \overrightarrow{\mathrm{ERP}}\right)+\mathrm{ZTD}_{\mathrm{r}} \cdot \mathrm{map}+ \\ & \mathrm{AMB}_{\mathrm{r}}^{\mathrm{s}}+\operatorname{clk}_{\mathrm{r}}\left(t_i\right)+\operatorname{clk}_{\mathrm{s}}\left(t_i\right)+ \\ & f\left(\vec{X}_{0 \mathrm{~s}}, \vec{X}_{0 \mathrm{r}}, \overrightarrow{\mathrm{ERP}}_0\right)+\varepsilon\end{aligned}$

(2)

式中，$\vec{X}_{0 \mathrm{~s}}、\vec{X}_{0 \mathrm{r}}、\overrightarrow{\mathrm{ERP}}_0$分别为$\vec{X}_{\mathrm{s}}、\vec{X}_{\mathrm{r}}、\overrightarrow{\mathrm{ERP}}$的近似值；$\Delta \vec{X}_{\mathrm{s}}、\Delta \vec{X}_{\mathrm{r}}、\Delta \overrightarrow{\mathrm{ERP}}$为待估的、相对于近似值的改正量。方程中待估参数的有效期见表 2，其中n_epoch为历元数，n_sta为地面站个数，n_sat为卫星数。

由表 2可以看出，卫星钟差参数和接收机钟差参数个数总量远大于其他待估参数个数。然而钟差参数的有效期仅为一个历元，属于局部参数，因此将其约化处理进行定轨解算。式(2)的向量形式可表示为：

$\begin{array}{c} P_j\left(t_i\right)= \\ {\left[\begin{array}{lllll}\frac{\partial f}{\partial \Delta \vec{X}_{\mathrm{s}}} & \frac{\partial f}{\partial \Delta \vec{X}_{\mathrm{r}}} & \frac{\partial f}{\partial \Delta \overrightarrow{\mathrm{ERP}}} & \text { map } & 1\end{array}\right]\left[\begin{array}{c}\Delta \vec{X}_{\mathrm{s}} \\ \Delta \vec{X}_{\mathrm{r}} \\ \Delta \overrightarrow{\mathrm{ERP}} \\ \mathrm{ZTD}_{\mathrm{r}} \\ \mathrm{AMB}_{\mathrm{r}}^{\mathrm{s}}\end{array}\right]+} \\ {\left[\begin{array}{ll}-1 & 1\end{array}\right]\left[\begin{array}{l}\mathrm{clk}_{\mathrm{s}}\left(t_i\right) \\ \operatorname{clk}_{\mathrm{r}}\left(t_i\right)\end{array}\right]+f\left(\vec{X}_{0 \mathrm{~s}}, \vec{X}_{0 \mathrm{r}}, \overrightarrow{\mathrm{ERP}}_0\right)+\varepsilon} \end{array}$

(3)

更一般地，历元i第j条观测的观测向量方程可表示为：

$\left\{ \begin{array}{l} {Y_{i, j}} = {\mathit{\boldsymbol{B}}_{i, j}}\mathit{\boldsymbol{X}} = \left[ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{B}}_{Gi, j}}}&{{\mathit{\boldsymbol{B}}_{Ri, j}}} \end{array}} \right]\left[ {{{\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{X}}_G}}&{\mathit{\boldsymbol{X}}_R^i} \end{array}}^{\rm{T}}}} \right.\\ {{\bf{B}}_{Gi, j}} = \left[ {\frac{{\partial f}}{{\partial \Delta {{\vec X}_{\rm{s}}}}}\;\;\; \cdots \;\;\;\frac{{\partial f}}{{\partial \Delta {{\vec X}_{\rm{r}}}}}\;\;\; \cdots } \right.\\ \left. {\;\;\;\;\;\frac{{\partial f}}{{\partial \Delta \overrightarrow {{\rm{ERP}}} }}\;\;\;{\mathop{\rm map}\nolimits} \;\;\;1\;\;\; \cdots } \right]\\ {\mathit{\boldsymbol{B}}_{Ri, j}} = \left[ {\begin{array}{*{20}{l}} \cdots &{ - 1}&1& \cdots \end{array}} \right]\\ {\mathit{\boldsymbol{X}}_G} = \left[ {\begin{array}{*{20}{l}} {\Delta {{\vec X}_{\rm{s}}}}& \cdots &{\Delta {{\vec X}_{\rm{r}}}}& \cdots &{\Delta \overrightarrow {{\rm{ERP}}} } \end{array}} \right.\\ \;\;\;{\left. {\begin{array}{*{20}{l}} {{\rm{ZT}}{{\rm{D}}_{\rm{r}}}}&{{\rm{AMB}}_{\rm{r}}^{\rm{s}}}& \cdots \end{array}} \right]^{\rm{T}}}\\ \mathit{\boldsymbol{X}}_R^i = {\left[ {\begin{array}{*{20}{l}} \cdots &{{{{\mathop{\rm clk}\nolimits} }_{\rm{s}}}\left( {{t_i}} \right)}&{{\rm{cl}}{{\rm{k}}_{\rm{r}}}\left( {{t_i}} \right)}& \cdots \end{array}} \right]^{\rm{T}}}\\ {Y_{i, j}} = {P_j}\left( {{t_i}} \right) - f\left( {{{\vec X}_{0{\rm{s}}}}, {{\vec X}_{0{\rm{r}}}}, {{\overrightarrow {{\rm{ERP}}} }_0}} \right) \end{array} \right.$

(4)

式中，系数矩阵B_{i, j}=[B_{Gi, j}  B_{Ri, j}]为历元i第j条观测的系数矩阵，其中B_{Gi, j}为该观测的全局变量系数矩阵，B_{Ri, j}为该观测的局部变量系数矩阵；待估参数向量X=[X_G  X_Rⁱ]^T，其中X_G为所有待估全局待变量，X_Rⁱ为仅在i历元时刻有效的待估参数；Y_{i, j}为观测常值。

基于式(4)、最小二乘原理及矩阵运算原理可知，历元i法方程N_bbⁱX=Wⁱ中，矩阵N_bbⁱ和Wⁱ可分别对历元i各观测法方程系数矩阵N_bb^ij和右矩阵W^ij叠加计算获得：

$\begin{gathered} \boldsymbol{N}_{b b}^i=\sum\limits_{j=1}^m \boldsymbol{N}_{b b}^{i j}=\sum\limits_{j=1}^m\left(\boldsymbol{B}_{i, j}^{\mathrm{T}} \boldsymbol{B}_{i, j}\right)= \\ \sum\limits_{j=1}^m\left[\begin{array}{ll} \boldsymbol{B}_{G i, j}^{\mathrm{T}} \boldsymbol{B}_{G i, j} & \boldsymbol{B}_{G i, j}^{\mathrm{T}} \boldsymbol{B}_{R i, j} \\ \boldsymbol{B}_{R i, j}^{\mathrm{T}} \boldsymbol{B}_{G i, j} & \boldsymbol{B}_{R i, j}^{\mathrm{T}} \boldsymbol{B}_{R i, j} \end{array}\right]=\\ {\left[\begin{array}{cc} \sum\limits_{j=1}^m\left(\boldsymbol{B}_{G i, j}^{\mathrm{T}} \boldsymbol{B}_{G i, j}\right) & \sum\limits_{j=1}^m\left(\boldsymbol{B}_{G i, j}^{\mathrm{T}} \boldsymbol{B}_{R i, j}\right) \\ \sum\limits_{j=1}^m\left(\boldsymbol{B}_{R i, j}^{\mathrm{T}} \boldsymbol{B}_{G i, j}\right) & \sum\limits_{j=1}^m\left(\boldsymbol{B}_{R i, j}^{\mathrm{T}} \boldsymbol{B}_{R i, j}\right) \end{array}\right]=} \\ {\left[\begin{array}{ll} \boldsymbol{N}_{b b G G i} & \boldsymbol{N}_{b b G R i} \\ \boldsymbol{N}_{b b R G i} & \boldsymbol{N}_{b b R R i} \end{array}\right]} \end{gathered}$

(5)

$\begin{gathered} \boldsymbol{W}^i=\sum\limits_{j=1}^m \boldsymbol{W}^{i j}=\sum\limits_{j=1}^m\left(\boldsymbol{B}_{i, j}^{\mathrm{T}} Y_{i, j}\right)= \\ \sum\limits_{j=1}^m\left[\begin{array}{l} \boldsymbol{B}_{G i, j}^{\mathrm{T}} Y_{i, j} \\ \boldsymbol{B}_{R i, j}^{\mathrm{T}} Y_{i, j} \end{array}\right]=\left[\begin{array}{l} \sum\limits_{j=1}^m \boldsymbol{B}_{G i, j}^{\mathrm{T}} Y_{i, j} \\ \sum\limits_{j=1}^m \boldsymbol{B}_{R i, j}^{\mathrm{T}} Y_{i, j} \end{array}\right]=\left[\begin{array}{l} \boldsymbol{W}_G^i \\ \boldsymbol{W}_R^i \end{array}\right] \end{gathered}$

(6)

式中，m为历元i中的观测数量。根据式(5)和(6)，历元i的法方程可写为：

$\left[\begin{array}{ll} \boldsymbol{N}_{b b G G i} & \boldsymbol{N}_{b b G R i} \\ \boldsymbol{N}_{b b R G i} & \boldsymbol{N}_{b b R R i} \end{array}\right]\left[\begin{array}{l} \boldsymbol{X}_G \\ \boldsymbol{X}_R^i \end{array}\right]=\left[\begin{array}{l} \boldsymbol{W}_G^i \\ \boldsymbol{W}_R^i \end{array}\right]$

(7)

根据式(7)可得：

$\boldsymbol{X}_R^i=\left(\boldsymbol{N}_{b b R R i}\right)^{-1}\left(\boldsymbol{W}_R^i-\boldsymbol{N}_{b b R G i} \boldsymbol{X}_G\right)$

(8)

将式(8)代入式(7)，可获得消去局部钟差参数后的历元i的法方程：

$\left\{\begin{array}{l} \boldsymbol{N}_{b b G G i^{\prime}} \boldsymbol{X}_G=\boldsymbol{W}_{G^{\prime}}^i \\ \boldsymbol{N}_{b b G G i^{\prime}}=\left(\boldsymbol{N}_{b b G G i}-\boldsymbol{N}_{b b G R i}\left(\boldsymbol{N}_{b b R R i}\right)^{-1} \boldsymbol{N}_{b b R G i}\right) \\ \boldsymbol{W}_{G^{\prime}}^i=\boldsymbol{W}_G^i-\boldsymbol{N}_{b b G R i}\left(\boldsymbol{N}_{b b R R i}\right)^{-1} \boldsymbol{W}_R^i \end{array}\right.$

(9)

根据式(9)对各历元法方程进行叠加，可建立所有历元观测数据统一的法方程：

$\left.\sum\limits_{i=1}^n \boldsymbol{N}_{b b G G i^{\prime}}\right) \boldsymbol{X}_G=\sum\limits_{i=1}^n \boldsymbol{W}_{G^{\prime}}^i$

(10)

式中，n为历元个数。对上式求逆解算可得：

$\left.\boldsymbol{X}_G=\sum\limits_{i=1}^n \boldsymbol{N}_{b b G G i^{\prime}}\right)^{-1} \quad \sum\limits_{i=1}^n \boldsymbol{W}_{G^{\prime}}^i$

(11)

式中，X_G为最小二乘全局参数最优解。将X_G代入式(8)可得到各历元时刻的局部参数解。

1.2 解算耗时分析

通过软件进行导航卫星精密定轨解算是基于式(5)~(11)分步实现的，其实质是对矩阵N_bbGGi、N_bbRRi、N_bbGRi、N_bbRGi、(N_bbRRi)^－1、N_bbGRi(N_bbRRi)^－1N_bbRGi和$\left(\sum\limits_{i=1}^n \boldsymbol{N}_{b b GG i^{\prime}}{ }^{-1}\right.$依次进行求解。表 3为这些矩阵运算的复杂度分析，其中n_g为全局待估参数个数，n_r为局部待估参数个数，n_epoch为历元数，m为一个历元中的观测数量。软件中采用Cholesky分解求逆。根据上节解算原理可知，表 3中序号1~6的矩阵每个历元都需要解算，所有历元的计算复杂度可表示为n_epoch×O(O表示单个历元矩阵运算的复杂度)。在大量采用地面站和卫星数的情况下，n_g≥n_r。从表 3可知，N_bbGRi(N_bbRRi)^－1N_bbRGi和$\left.\sum\limits_{i=1}^n \boldsymbol{N}_{b b G Gi^{\prime}}\right)^{-1}$的计算最耗时，分别为式(9)中的钟差参数约化和式(11)中的统一法方程矩阵求逆。因此，本文仅针对二者进行优化，以提高导航卫星精密定轨解算效率。

2 GNSS卫星精密定轨解算效率优化

为充分发挥ARM架构CPU的计算效能，采用基于CPU的多核、多线程并行运算技术和基于CPU的分层存储架构对上述2个计算耗时最长的矩阵进行优化。采用多线程方法对钟差参数约化解算部分进行优化，采用OpenBlas对统一方法求逆部分进行优化。OpenBlas是同时采用上述2种方法实现矩阵运算效率提升的开源高性能科学计算库。

2.1 基于多线程的消参并行解算优化

通过式(5)可以获得N_bbGRi、N_bbGRi和N_bbRRi，然后计算N_bbGRi(N_bbRRi)^－1N_bbRGi，最后完成式(9)的计算。由此可知，钟差参数约化解算部分N_bbGRi(N_bbRRi)^－1N_bbRGi在解算过程中是完全独立的，可以单独对其进行优化。优化过程为：首先优化矩阵运算N_bbGRi(N_bbRRi)^－1，然后优化矩阵运算N_bbGRi(N_bbRRi)^－1N_bbRGi，这2步的优化方法一致。矩阵乘法运算可表示为C_I×J=A_I×PB_P×J。软件中代码实现是对浮点数的一个3层循环运算，见图 1，图中a_ip、b_pj、c_ij分别为矩阵A、B、C的元素，是一个串行执行过程。

根据矩阵乘法特点，可对图 1中最外2层循环(即矩阵A的行和矩阵B的列)进行拆分，分解为多个独立的运算任务，进而采用多线程并行优化。假如采用n(n≤I)个线程进行并行运算，对每个线程分配最外层循环个数可表示为：$\left\{\begin{array}{l}I_l=S=\operatorname{ceil}(I / n) \\ I_n=I-(n-1) S\end{array}(l=1, 2, \cdots, n-1\right.$; ceil表示向上取整)。按照最理想的情况计算，并行后运算效率为原来单线程的n倍。但实际上，n取值需要权衡计算机CPU资源，n越大，分配给每个线程的CPU缓存资源越少，可能会增加数据在存储单元上的传输时间。假如在n个线程情况下数据的传输时间为原来单线程的m倍，则n个线程实际的运行效率为原来单线程的n/m倍。因此，需要权衡计算机CPU资源和运算数据量，使n/m值最大，以获取最高的并行解算效率。

2.2 基于OpenBlas的矩阵求逆解算优化

通过式(9)和(10)分别获得矩阵$\sum\limits_{i=1}^n \boldsymbol{N}_{b b G G i^{\prime}}$、$\sum\limits_{i=1}^n \boldsymbol{W}_{G^{\prime}}^i$，然后计算$\left(\sum\limits_{i=1}^n \boldsymbol{N}_{b b G Gi^{\prime}}\right)^{-1}$，最后完成式(11)的解算。由此可知，统一法方程求逆部分$\left.\sum\limits_{i=1}^n \boldsymbol{N}_{b b G G i^{\prime}}\right)^{-1}$在解算过程中是完全独立的，同样可以单独对其进行解算优化。软件中对$\left.\sum\limits_{i=1}^n \boldsymbol{N}_{b b G G i^{\prime}}\right)^{-1}$的解算采用Cholesky分解，其中矩阵运算采用Lapack库，而本文采用OpenBlas库。

Lapack相比基础线性代数程序从算法层面进行了许多优化，OpenBlas继承了这些优点，并进一步采用前文基于CPU技术特点的2个优化方法对矩阵运算进行优化。基于CPU分层存储架构优化计算效率的具体思路见文献[9]。

在多线程并行优化方面，OpenBlas采用仅将矩阵B作列拆分的方式进行多线程优化^[10]。在此基础上，考虑每个线程独有的缓存资源以及共享缓存资源，根据线程数调整内核GEBP^[9]的大小，使内核GEBP适合缓存大小，以减少缓存不命中的情况，详细过程见文献[10]。

3 实验分析

CPU具体架构信息见表 4，GNSS卫星精密定轨软件的具体解算策略见表 5。采用32颗卫星和不同数目的测站检验精密定轨解算优化情况，每次定轨解算执行8次迭代，每次迭代均重复执行消参和法方程求逆解算。

在飞腾CPU上进行解算时，采用16个线程的多线程方法对消参部分进行优化，优化前后解算效率对比见图 2。可以看出，相比于原始单线程，采用多线程后消参解算效率大幅提升。随着测站数量的增加，采用单线程消参解算耗时急剧上升，而多线程解算耗时平缓上升，且多线程解算效率提升效果更明显。在100个测站的情况下，统计第1次迭代中各历元单线程和多线程消参解算耗时。单线程解算1个历元平均耗时1.105 s，而多线程仅为0.188 s。100个测站情况下所有8次迭代消参耗时情况见图 3。由图可知，与第1次迭代耗时情况一致，多线程解算效率约为单线程解算效率的6倍。

对于求逆解算，将采用单线程和多线程矩阵运算优化方法的OpenBlas库和原始Lapack库进行解算对比。OpenBlas库采用16线程，统计每次解算中8次迭代统一法方程求逆耗时的平均值，图 4为不同测站情况下采用Lapack库和OpenBlas库法方程求逆平均耗时对比。可以看出，OpenBlas库解算效率明显优于Lapack库，且随着测站数目的增加，OpenBlas库的优势更加明显。图 5为在100个测站情况下，8次迭代中二者法方程求逆解算耗时对比。经统计，采用Lapack库平均耗时为2 264 s，采用OpenBlas库平均耗时仅为78 s。

表 6为本文优化策略下使用ARM CPU和Intel CPU单次定轨的平均迭代时间。由表可知，在ARM CPU上通过上述2种方法优化(ARM+多线程+OpenBlas)的定轨解算效率为表 1中原始定轨处理(ARM+Lapack)的6倍左右，并且超过原始软件在Intel CPU(Intel+MKL)上的解算效率。本文优化方法在Intel CPU上也可以提高解算效率，采用MKL(Intel+多线程+MKL)时解算效率还是最高的。ARM CPU上最优解算耗时与Intel CPU上最优解算耗时相比，从原来的5倍缩小到3倍，说明本文方法可以提升ARM CPU的计算效能，提高导航卫星定轨解算效率。

4 结语

以国产飞腾CPU为例，讨论在国产ARM架构CPU基础上的导航卫星精密定轨解算效率优化方法。影响导航卫星精密定轨解算效率的最主要因素为钟差参数约化和法方程求逆运算，对前者采用多线程优化方法，对后者采用OpenBlas开源库进行优化。优化后，2个部分的解算效率均大幅提升，且随着所需解算参数的增加，效率提升的倍数不断增加。在100个测站32颗卫星和OpenBlas采用16线程的情况下，针对矩阵乘法同时进行单线程和多线程运算优化的法方程求逆解算时，OpenBlas库的解算效率为Lapack库的29倍，这个优化倍数已突破仅通过多线程优化所能达到的最理想倍数(16倍)，说明OpenBlas中单线程矩阵运算优化同样可以大幅提升解算效率。后续研究将进一步根据ARM CPU的架构特点和OpenBlas的原理，针对ARM架构对OpenBlas内部参数进行适配性调整，以进一步提升ARM CPU解算卫星精密轨道的效能。