基于國產ARM架構CPU的導航衛星精密定軌解算效率優化方法

廖 敏 唐成盼 周善石 陳建兵 胡小工 馮學斌 陳桂根 李 凱

1 中國科學院上海天文臺,上海市南丹路80號,200030

2 中國科學院大學,北京市玉泉路19號甲,100049

3 探索數據科技(深圳)有限公司,深圳市留仙洞路33號,518055

提高導航衛星精密定軌解算效率是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還需要尋求更合適的優化方法以提升導航衛星精密定軌的解算效率。

表1 Intel CPU和ARM CPU單次迭代平均耗時

本文選取國產飛騰CPU的最新型號FT-2000+進行衛星精密定軌解算效率優化研究。

1 GNSS衛星精密定軌解算

1.1 解算原理

設歷元i的第j條觀測是接收機r與衛星s的偽距相位觀測,添加動力學約束的觀測方程為:

(1)

將式(1)在近似值處線性化得到:

(2)

表2 GNSS精密定軌中待估參數的時效性

由表2可以看出,衛星鐘差參數和接收機鐘差參數個數總量遠大于其他待估參數個數。然而鐘差參數的有效期僅為一個歷元,屬于局部參數,因此將其約化處理進行定軌解算。式(2)的向量形式可表示為:

Pj(ti)=

(3)

更一般地,歷元i第j條觀測的觀測向量方程可表示為:

(4)

(5)

(6)

式中,m為歷元i中的觀測數量。根據式(5)和(6),歷元i的法方程可寫為:

(7)

根據式(7)可得:

(8)

將式(8)代入式(7),可獲得消去局部鐘差參數后的歷元i的法方程:

(9)

根據式(9)對各歷元法方程進行疊加,可建立所有歷元觀測數據統一的法方程:

(10)

式中,n為歷元個數。對上式求逆解算可得:

(11)

式中,XG為最小二乘全局參數最優解。將XG代入式(8)可得到各歷元時刻的局部參數解。

1.2 解算耗時分析

表3 導航衛星精密定軌解算各環節復雜度分析

2 GNSS衛星精密定軌解算效率優化

為充分發揮ARM架構CPU的計算效能,采用基于CPU的多核、多線程并行運算技術和基于CPU的分層存儲架構對上述2個計算耗時最長的矩陣進行優化。采用多線程方法對鐘差參數約化解算部分進行優化,采用OpenBlas對統一方法求逆部分進行優化。OpenBlas是同時采用上述2種方法實現矩陣運算效率提升的開源高性能科學計算庫。

2.1 基于多線程的消參并行解算優化

通過式(5)可以獲得NbbGRi、NbbGRi和NbbRRi,然后計算NbbGRi(NbbRRi)-1NbbRGi,最后完成式(9)的計算。由此可知,鐘差參數約化解算部分NbbGRi(NbbRRi)-1NbbRGi在解算過程中是完全獨立的,可以單獨對其進行優化。優化過程為:首先優化矩陣運算NbbGRi(NbbRRi)-1,然后優化矩陣運算NbbGRi(NbbRRi)-1NbbRGi,這2步的優化方法一致。矩陣乘法運算可表示為CI×J=AI×PBP×J。軟件中代碼實現是對浮點數的一個3層循環運算,見圖1,圖中aip、bpj、cij分別為矩陣A、B、C的元素,是一個串行執行過程。

圖1 矩陣乘法的3層循環實現

2.2 基于OpenBlas的矩陣求逆解算優化

Lapack相比基礎線性代數程序從算法層面進行了許多優化,OpenBlas繼承了這些優點,并進一步采用前文基于CPU技術特點的2個優化方法對矩陣運算進行優化。基于CPU分層存儲架構優化計算效率的具體思路見文獻[9]。

在多線程并行優化方面,OpenBlas采用僅將矩陣B作列拆分的方式進行多線程優化[10]。在此基礎上,考慮每個線程獨有的緩存資源以及共享緩存資源,根據線程數調整內核GEBP[9]的大小,使內核GEBP適合緩存大小,以減少緩存不命中的情況,詳細過程見文獻[10]。

3 實驗分析

CPU具體架構信息見表4,GNSS衛星精密定軌軟件的具體解算策略見表5。采用32顆衛星和不同數目的測站檢驗精密定軌解算優化情況,每次定軌解算執行8次迭代,每次迭代均重復執行消參和法方程求逆解算。

表4 國產飛騰CPU和Intel CPU架構信息

表5 GNSS衛星精密定軌軟件解算策略

在飛騰CPU上進行解算時,采用16個線程的多線程方法對消參部分進行優化,優化前后解算效率對比見圖2。可以看出,相比于原始單線程,采用多線程后消參解算效率大幅提升。隨著測站數量的增加,采用單線程消參解算耗時急劇上升,而多線程解算耗時平緩上升,且多線程解算效率提升效果更明顯。在100個測站的情況下,統計第1次迭代中各歷元單線程和多線程消參解算耗時。單線程解算1個歷元平均耗時1.105 s,而多線程僅為0.188 s。100個測站情況下所有8次迭代消參耗時情況見圖3。由圖可知,與第1次迭代耗時情況一致,多線程解算效率約為單線程解算效率的6倍。

圖2 消參部分采用單線程和多線程解算耗時對比

圖3 100個測站解算8次迭代消參耗時

對于求逆解算,將采用單線程和多線程矩陣運算優化方法的OpenBlas庫和原始Lapack庫進行解算對比。OpenBlas庫采用16線程,統計每次解算中8次迭代統一法方程求逆耗時的平均值,圖4為不同測站情況下采用Lapack庫和OpenBlas庫法方程求逆平均耗時對比。可以看出,OpenBlas庫解算效率明顯優于Lapack庫,且隨著測站數目的增加,OpenBlas庫的優勢更加明顯。圖5為在100個測站情況下,8次迭代中二者法方程求逆解算耗時對比。經統計,采用Lapack庫平均耗時為2 264 s,采用OpenBlas庫平均耗時僅為78 s。

圖4 采用Lapack庫和OpenBlas庫法方程求逆平均耗時對比

圖5 8次迭代法方程求逆耗時對比

表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的計算效能,提高導航衛星定軌解算效率。

表6 優化后ARM CPU和Intel CPU 單次定軌的平均迭代時間

4 結 語

以國產飛騰CPU為例,討論在國產ARM架構CPU基礎上的導航衛星精密定軌解算效率優化方法。影響導航衛星精密定軌解算效率的最主要因素為鐘差參數約化和法方程求逆運算,對前者采用多線程優化方法,對后者采用OpenBlas開源庫進行優化。優化后,2個部分的解算效率均大幅提升,且隨著所需解算參數的增加,效率提升的倍數不斷增加。在100個測站32顆衛星和OpenBlas采用16線程的情況下,針對矩陣乘法同時進行單線程和多線程運算優化的法方程求逆解算時,OpenBlas庫的解算效率為Lapack庫的29倍,這個優化倍數已突破僅通過多線程優化所能達到的最理想倍數(16倍),說明OpenBlas中單線程矩陣運算優化同樣可以大幅提升解算效率。后續研究將進一步根據ARM CPU的架構特點和OpenBlas的原理,針對ARM架構對OpenBlas內部參數進行適配性調整,以進一步提升ARM CPU解算衛星精密軌道的效能。