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

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

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

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

提高導(dǎo)航衛(wèi)星精密定軌解算效率是GNSS衛(wèi)星精密定軌技術(shù)的研究重點(diǎn)之一[1-2]。Ge等[3]采用消參方法,在每個(gè)歷元時(shí)刻消去法方程中的過時(shí)參數(shù),極大地減少了最終法方程解算參數(shù)的個(gè)數(shù),從而有效提高了導(dǎo)航衛(wèi)星精密定軌解算的效率,該方法是從算法層面對解算效率進(jìn)行優(yōu)化。隨著CPU的迭代發(fā)展,部分學(xué)者開展基于CPU技術(shù)提高計(jì)算效率的研究,包括基于CPU的多核、多線程并行運(yùn)算技術(shù)[4-5]和基于CPU分層存儲架構(gòu)的方法[6]。

目前,關(guān)于提高GNSS數(shù)據(jù)處理軟件計(jì)算效率的研究幾乎都是基于X86架構(gòu)的Intel CPU,涉及國產(chǎn)ARM架構(gòu)CPU[7]上效率提升的實(shí)驗(yàn)和研究很少。在科學(xué)計(jì)算中,線性代數(shù)庫是運(yùn)行效率優(yōu)化的重要工具,Lapack能以高效的方式實(shí)現(xiàn)科學(xué)計(jì)算的接口規(guī)范[8],并且適用于X86、ARM和PowerPC等不同架構(gòu)平臺。Intel CPU發(fā)展至今已經(jīng)很成熟,Intel可實(shí)現(xiàn)適應(yīng)其CPU的高效運(yùn)算庫Intel Math Kernel Library(MKL)。表1為導(dǎo)航衛(wèi)星精密定軌處理軟件分別運(yùn)行于Intel CPU中性能突出的X86架構(gòu)Intel(R) Xeon(R) Gold 6134 CPU(表中Intel)和國產(chǎn)CPU中性能突出的飛騰ARM FT-2000+CPU(表中ARM)的運(yùn)行耗時(shí),軟件運(yùn)行時(shí)采用相同的數(shù)據(jù)和場景?？梢钥闯?采用Lapack時(shí),Intel CPU的解算效率是ARM CPU的3倍多;采用MKL后,Intel CPU的解算效率是ARM CPU的5倍多。由此可知,ARM CPU還需要尋求更合適的優(yōu)化方法以提升導(dǎo)航衛(wèi)星精密定軌的解算效率。

表1 Intel CPU和ARM CPU單次迭代平均耗時(shí)

本文選取國產(chǎn)飛騰CPU的最新型號FT-2000+進(jìn)行衛(wèi)星精密定軌解算效率優(yōu)化研究。

1 GNSS衛(wèi)星精密定軌解算

1.1 解算原理

設(shè)歷元i的第j條觀測是接收機(jī)r與衛(wèi)星s的偽距相位觀測,添加動(dòng)力學(xué)約束的觀測方程為:

(1)

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

(2)

表2 GNSS精密定軌中待估參數(shù)的時(shí)效性

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

Pj(ti)=

(3)

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

(4)

(5)

(6)

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

(7)

根據(jù)式(7)可得:

(8)

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

(9)

根據(jù)式(9)對各歷元法方程進(jìn)行疊加,可建立所有歷元觀測數(shù)據(jù)統(tǒng)一的法方程:

(10)

式中,n為歷元個(gè)數(shù)。對上式求逆解算可得:

(11)

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

1.2 解算耗時(shí)分析

表3 導(dǎo)航衛(wèi)星精密定軌解算各環(huán)節(jié)復(fù)雜度分析

2 GNSS衛(wèi)星精密定軌解算效率優(yōu)化

為充分發(fā)揮ARM架構(gòu)CPU的計(jì)算效能,采用基于CPU的多核、多線程并行運(yùn)算技術(shù)和基于CPU的分層存儲架構(gòu)對上述2個(gè)計(jì)算耗時(shí)最長的矩陣進(jìn)行優(yōu)化。采用多線程方法對鐘差參數(shù)約化解算部分進(jìn)行優(yōu)化,采用OpenBlas對統(tǒng)一方法求逆部分進(jìn)行優(yōu)化。OpenBlas是同時(shí)采用上述2種方法實(shí)現(xiàn)矩陣運(yùn)算效率提升的開源高性能科學(xué)計(jì)算庫。

2.1 基于多線程的消參并行解算優(yōu)化

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

圖1 矩陣乘法的3層循環(huán)實(shí)現(xiàn)

2.2 基于OpenBlas的矩陣求逆解算優(yōu)化

Lapack相比基礎(chǔ)線性代數(shù)程序從算法層面進(jìn)行了許多優(yōu)化,OpenBlas繼承了這些優(yōu)點(diǎn),并進(jìn)一步采用前文基于CPU技術(shù)特點(diǎn)的2個(gè)優(yōu)化方法對矩陣運(yùn)算進(jìn)行優(yōu)化。基于CPU分層存儲架構(gòu)優(yōu)化計(jì)算效率的具體思路見文獻(xiàn)[9]。

在多線程并行優(yōu)化方面,OpenBlas采用僅將矩陣B作列拆分的方式進(jìn)行多線程優(yōu)化[10]。在此基礎(chǔ)上,考慮每個(gè)線程獨(dú)有的緩存資源以及共享緩存資源,根據(jù)線程數(shù)調(diào)整內(nèi)核GEBP[9]的大小,使內(nèi)核GEBP適合緩存大小,以減少緩存不命中的情況,詳細(xì)過程見文獻(xiàn)[10]。

3 實(shí)驗(yàn)分析

CPU具體架構(gòu)信息見表4,GNSS衛(wèi)星精密定軌軟件的具體解算策略見表5。采用32顆衛(wèi)星和不同數(shù)目的測站檢驗(yàn)精密定軌解算優(yōu)化情況,每次定軌解算執(zhí)行8次迭代,每次迭代均重復(fù)執(zhí)行消參和法方程求逆解算。

表4 國產(chǎn)飛騰CPU和Intel CPU架構(gòu)信息

表5 GNSS衛(wèi)星精密定軌軟件解算策略

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

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

圖3 100個(gè)測站解算8次迭代消參耗時(shí)

對于求逆解算,將采用單線程和多線程矩陣運(yùn)算優(yōu)化方法的OpenBlas庫和原始Lapack庫進(jìn)行解算對比。OpenBlas庫采用16線程,統(tǒng)計(jì)每次解算中8次迭代統(tǒng)一法方程求逆耗時(shí)的平均值,圖4為不同測站情況下采用Lapack庫和OpenBlas庫法方程求逆平均耗時(shí)對比?？梢钥闯?OpenBlas庫解算效率明顯優(yōu)于Lapack庫,且隨著測站數(shù)目的增加,OpenBlas庫的優(yōu)勢更加明顯。圖5為在100個(gè)測站情況下,8次迭代中二者法方程求逆解算耗時(shí)對比。經(jīng)統(tǒng)計(jì),采用Lapack庫平均耗時(shí)為2 264 s,采用OpenBlas庫平均耗時(shí)僅為78 s。

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

圖5 8次迭代法方程求逆耗時(shí)對比

表6為本文優(yōu)化策略下使用ARM CPU和Intel CPU單次定軌的平均迭代時(shí)間。由表可知,在ARM CPU上通過上述2種方法優(yōu)化(ARM+多線程+OpenBlas)的定軌解算效率為表1中原始定軌處理(ARM+Lapack)的6倍左右,并且超過原始軟件在Intel CPU(Intel+MKL)上的解算效率。本文優(yōu)化方法在Intel CPU上也可以提高解算效率,采用MKL(Intel+多線程+MKL)時(shí)解算效率還是最高的。ARM CPU上最優(yōu)解算耗時(shí)與Intel CPU上最優(yōu)解算耗時(shí)相比,從原來的5倍縮小到3倍,說明本文方法可以提升ARM CPU的計(jì)算效能,提高導(dǎo)航衛(wèi)星定軌解算效率。

表6 優(yōu)化后ARM CPU和Intel CPU 單次定軌的平均迭代時(shí)間

4 結(jié) 語

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