高軌導航接收機的時鐘輔助定位算法

郭坤昊, 葛 建, 陳夏蘭, 房志博, 高 銘

(1. 中國科學院空天信息創新研究院, 北京 100094;2. 中國科學院大學電子電氣與通信工程學院, 北京 100049)

0 引言

運行在地球靜止軌道(geostationary orbit satellite, GEO)、傾斜地球同步軌道 (inclined geo-synchronous orbit,IGSO)的衛星的定軌通常需要依靠地面的測控網絡實現。這種定軌方式成本高并且時效性差,如果能像運行在低地球軌道(low-Earth orbit,LEO)上的衛星那樣,通過在星上搭載全球衛星導航系統(global navigation satellite system,GNSS)接收機,實現對衛星的定軌,將能為衛星發射和運行帶來極大的便利[1]。適合在GEO、IGSO等衛星上工作的GNSS接收機屬于高軌GNSS接收機,其典型工作高度高于一般導航衛星的軌道高度。

高軌GNSS接收機接收的導航衛星信號條件與地面和LEO衛星上的信號條件完全不同,主要表現在信號弱和幾何構形差兩個方面。

首先,GNSS衛星播發的導航信號是朝向地球發射的,全球定位系統(global positioning system,GPS)L1信號的雙邊主波束寬度大約是32°,僅能夠覆蓋地球及1000 km以下LEO衛星軌道[2]。高軌GNSS接收機恰好位于主瓣照射區域的機會非常小(如圖1的紅色區域)。如果僅使用導航衛星的主波束進行定位,即使用上GNSS星座的所有導航衛星,可用衛星數也很難達到4顆。為了有更多可用的導航衛星,高軌GNSS接收機的靈敏度通常做得很高,從而可以使用部分導航衛星的第1、第2旁瓣進行定位。但是,導航衛星的第1、第2旁瓣的信號功率在三維天線方向圖的不同剖面上差異較大。即使在較強的剖面上,信號功率也要比主瓣低10 dB到20 dB[3]。此外,高軌接收機到導航衛星的距離大約是地面接收機的3倍左右,相應的自由空間傳播損耗增加約10 dB。在距離因素和導航衛星的天線方向圖因素的共同作用下,到達高軌GNSS接收機的信號功率比地面接收機低10 dB到30 dB[4]。

圖1 高軌GNSS接收機的接收信號條件Fig.1 Reception conditions of high orbit GNSS receivers

其次,由于導航衛星旁瓣信號的覆蓋區域,和地心形成的夾角最遠有60°,在位于GEO的接收機的視角上,所有可用的導航衛星都分布在地球兩側距地心不超過32°的狹小區域內(如圖1中高軌GNSS接收機與GNSS #2之間的相對關系),并且地心兩側10°的范圍內,因地球遮擋和電離層的影響沒有可用的導航衛星,這些導致接收機可見星的幾何構形較差。姜慧等人的仿真實驗表明對于運行在GEO軌道的高軌接收機利用北斗衛星導航系統(BeiDou satellite navigation system, BDS)定位時幾何精度因子值(geometric dilution of precision, GDOP)均值約為49,其最大值約161;位置精度因子值(po-sition dilution of precision,PDOP)的均值約為35,其最大值約115[5]。

上述高軌GNSS接收機的信號條件最終反映在接收機性能上的表現是定位精度低、連續性差。

由于信號弱,高軌GNSS接收機典型的接收信號載噪比C/N0值在30(dB·Hz)左右,對應的碼偽距測量精度大約在10 m左右,按PDOP值為35計算,由此導致的單歷元解算誤差大約350 m,遠大于地面接收機的誤差。同樣由于信號弱,能夠超過高軌GNSS接收機解調門限的導航衛星數較少。即使靈敏度較高的接收機,仍然存在可用衛星不足4顆的情況。秦紅磊等人的仿真結果表明,接收機靈敏度閾值設為-175 dBW時GEO衛星可見GPS衛星數目在99%的運行時間里都不足4顆;而當接收功率閾值降低至-190 dBW時,雖然可見星數目明顯增加,但2%的時間里可見星數目不足4顆[6]。

在較難改變信號接收功率和衛星幾何構形的情況下,為了解決高軌GNSS接收機定位精度低和連續性差的問題,一些學者提出了軌道外推、卡爾曼濾波和深組合導航等適用于高軌GNSS接收機的算法。學者詹鵬宇對高軌場景下的導航定軌采用了無味卡爾曼濾波算法,還結合龍格-庫塔方法的思想對擴展卡爾曼濾波算法進行了改進,并且對多種非線性濾波算法的高軌衛星定軌性能進行了對比分析[7];學者荊帥提出一種采用軌道外推器(orbit polating,OP)與GNSS觀測量進行深組合的方案,在測距域內重新構造GNSS/OP深組合濾波器,把傳統的解析式外推器或數值積分式外推器改進成一種簡化的基于經驗加速度模型輔助的四階龍格庫塔積分器,并選擇了GNSS空間服務域(sp-ace service volume,SSV)用戶-高橢圓軌道(highly elliptical orbit,HEO)飛行器作為測試對象,仿真驗證了GNSS/OP深組合方案的性能[8];學者林魁提出了融合 GNSS 與相對測量的序貫博弈自主定軌方法。基于SSV多星任務需求,利用高精度星間測量進一步提高航天器自主定軌精度,針對分布式系統中因狀態估計間的互相關關系導致濾波發散問題,提出了一種序貫博弈優化算法,并以四顆GEO衛星任務為例對算法進行了驗證[9]。

以上算法都存在運算量大的問題并且存在一定的適用條件。例如,學者詹鵬提出的利用卡爾曼濾波器作為時鐘“馴服-守時”模型的算法,雖然精度較高,但是算法實現時濾波器的協方差矩陣和狀態轉移矩陣等初始化存在困難;學者林魁提出的序貫博弈自主定軌方法由于需要星間鏈路支撐和多用戶之間建立星間鏈路,適用于多個用戶間高精度的定軌需求。本文提出利用接收機恒溫晶振的頻率穩定性抑制觀測噪聲,以較小的運算量顯著地改善高軌接收機定位的連續性和精度,圖2為流程圖。

圖2 算法流程圖Fig.2 Flowchart of the algorithm

1 時鐘輔助定位算法流程

1.1 算法概述

接收機硬件各通道跟蹤導航衛星信號所獲得的星上時觀測量ts(上標S表示導航衛星的序號)作為本算法的原始輸入,經圖2所示的流程計算,最終得到接收機位置矢量[x,y,z]的解算結果。圖2中“偽距構建”、“鐘差解算”和“時間維持”3個模塊構成了時間維持環路,用于產生歷元時間的優化預測值tk(下標k為觀測歷元的序號)。本算法需要用恒溫晶振(oven controlled crystal oscillator,OCXO)等穩定性高的頻率源維持歷元時間間隔的穩定。

1.2 偽距構建

接收機通過對測距碼相位的觀測獲取衛星信號離開導航衛星的時間,即星上時ts。接收機在第k個歷元同時觀測到N個導航衛星信號的發射時間構成觀測值輸入矢量[t1,t2,…,tN]k。根據接收機的歷元時刻tk和觀測到的星上時ts可構造第S顆衛星的偽距ρs,其表達式如下

ρs=c·(tk-ts)

(1)

其中,c為光速,所觀測的N顆衛星的偽距構建成偽距觀測值矢量[ρ1,ρ2,…,ρN]k。

1.3 鐘差解算

“鐘差解算模塊”用來估計tk與導航系統真實時間之間的偏差δtk。

在N≥4,并且時間精度因子值(time dilution of precision,TDOP)小于某個預先設定的閾值(比如60)的條件下,通過求解方程組(2),獲得δtk的估計值。

(2)

方程組(2)中,[xs,ys,zs]表示第S顆衛星的位置矢量,該矢量可以根據衛星星歷和信號發射時刻ts計算得到。[x,y,z]是待求的接收機位置矢量。由方程組(2)解出的接收機位置包含較多噪聲,當時間維持模塊進入穩態后,就不再使用方程組(2)解出的位置。

1.4 時間維持

1.4.1 時間維持模塊的正常運行狀態

當鐘差解算模塊輸出的δtk有效時,時間維持模塊需根據帶有較大噪聲的鐘差估計值δtk不斷修正本地鐘差和鐘漂參數,盡可能維持本地時間的準確[10]。

盡管可以采用高階比例積分微分控制器(proportion-integral-derivative,PID)算法[11]或者卡爾曼濾波算法[12]進一步提升性能,但實驗表明,對于充分老化的OCXO而言,簡單的二階鎖相環算法就能獲得相當滿意的效果[13]。具體迭代步驟如式(3)和式(4)所示

tk+1=tk+Ts+dk+C1·δtk

(3)

dk+1=dk+C2·δtk

(4)

其中,Ts代表觀測歷元時間間隔的標稱值,dk是第k歷元接收機鐘漂的估計值。

式(3)和式(4)中,參數C1和C2是二階鎖相環的環路參數[14],它們是由環路的阻尼系數和噪聲帶寬共同決定的。其中,C1和C2與環路自然頻率ωn、阻尼系數ξ的關系分別如式(5)和式(6)所示

(5)

(6)

而ωn和ξ又決定了環路的噪聲帶寬BL,其關系如式(7)

(7)

ξ通常取值為0.707,BL的選擇對算法性能的影響最為明顯,與GPS馴服恒溫晶振(GPSDO)類似,BL可以設計得非常小。考察方程組(2),由于衛星的高動態,等號左邊[xs,ys,zs]和[x,y,z]的真值都存在較大動態,而δtk的真值則僅包含由恒溫晶振阿倫方差導致的極小動態。等號右邊的偽距ρs則包含了衛星動態和恒溫晶振動態之和。接收機的碼跟蹤環需要維持對ρs的跟蹤,因此碼跟蹤環的噪聲帶寬Bn需要設置在2 Hz左右。經方程組(2)解出的δtk已剝離了衛星的高動態,所以時間維持環的噪聲帶寬可以做得很窄。利用GNSS模擬器模擬多種軌道場景的測試表明,如果采用OCXO做頻率源,即使將BL設為0.001 Hz,本算法仍可以很穩定地工作。

1.4.2 時間維持環路的初始化

式(3)和式(4)為迭代公式,初值d0可設為0,也可以將接收機上一次運行時dk的收斂值作為本次運行的初值。初值t0需要利用觀測到的某個導航星的星上時ts加上一個預估的傳播時延作為初始嘗試值,并利用第1次解得的鐘差值δt0修正得到t0的迭代初值。

1.4.3 時間維持模塊的守時狀態

迭代過程中,如果某個歷元可觀測衛星數N<4,或者TDOP大于設定的閾值,則鐘差觀測值δtk不可用,時間維持模塊按δtk=0繼續迭代運行。此時時間維持模塊處于守時狀態。在環路已收斂的條件下,OCXO的穩定性越好,則時間維持模塊可連續守時的時間就越長。

1.5 定位解算

時間維持模塊進入穩定狀態后,其輸出的歷元時間tk是對當前歷元的GNSS時的優化估計,通常比單點解算得到的GNSS時間更準確(見本文2.2節)。因此,按式(1)構造的偽距ρs非常接近接收機和導航衛星之間的真距。因此,在定位解算模塊中,通過求解不帶鐘差的方程組(8),得到關于接收機位置矢量[x,y,z]的優化估計。

(8)

應注意到,即使在N=3的情況下,方程組(8)依然能夠維持定位解算。

2 時鐘輔助算法的精度

2.1 高軌GNSS接收機的主要誤差源

GNSS接收機的定位誤差由測距誤差(user equivalent range error,UERE)和所觀測衛星的幾何構形共同決定。其中UERE又分為用戶測距誤差(user range error,URE)和用戶設備誤差(user equipment error,UEE)兩部分[16]。

URE主要包含導航衛星的星鐘誤差和軌道誤差以及電離層誤差和對流層誤差。對于高軌GNSS接收機而言,由于信號穿越電離層和對流層的時段占整個可視時段的比例很小,可以在選星策略中排除這些時段的觀測量,因此可以只考慮導航衛星的星鐘誤差和軌道誤差。這些誤差由衛星導航系統運控的精度決定。許多的實測數據表明鐘差和軌道誤差導致的URE精度在1.0 m(1σ)左右[17]。

URE的值是緩變的,在幾分鐘的時間范圍內,同一顆導航衛星的URE值幾乎為一個常數。

UEE主要包含接收機熱噪聲和多徑干擾。對于高軌 GNSS 接收機而言, 可以通過合理的天線布局等方法極大地減輕多徑的影響。因此,這里只考慮熱噪聲的影響。由熱噪聲導致的碼偽距測量誤差可由式(9)估算[18]。熱噪聲標準差σDLL和輸入信號載噪比的數學關系如式(9)所示

σDLL=λL·

(9)

其中,λL(m)表示測距碼碼片寬度;D(碼片數)表示相關器的間距;Bfe(Hz)表示雙邊前端帶寬;Tc表示碼片周期且Tc=1/Rc,Rc是碼片速率。由于高軌GNSS接收機接收的導航信號十分微弱,導致熱噪聲誤差可達10m左右。在高軌條件下,熱噪聲誤差是接收機最主要的誤差源。

UEE的值是快變的。如果接收機的歷元間隔大于1/(2·Bn),可以認為不同歷元間的熱噪聲測量誤差是相互獨立的,并且UEE與URE之間是相互獨立的。

此外,由熱噪聲導致的測量誤差UEE是0均值的。

2.2 tk誤差

與方程組(2)相對應,第k個歷元各星URE和UEE的瞬時值分別可構造成誤差列矢量UREk和UEEk,則第k個歷元的解算誤差[exk,eyk,ezk,C·eδtk]T可表示為式(10)

(10)

其中,G′為方程組(2)線性化后的系數矩陣。

方程(10)為線性方程,所以δtk的解算誤差序列eδtk可表示為由UREk導致的誤差序列eδtURE,k和由UEEk導致的誤差序列eδtUEE,k之和,如式(11)所示

eδtk=eδtURE,k+eδtUEE,k

(11)

由于各星URE值是緩變的,G′由衛星的幾何構形決定,也是緩變量,所以eδtURE,k是緩變量,其值在數分鐘內保持恒定。而UEE是0均值快變量,所以eδtUEE,k也是0均值的快變量。即eδtUEE,k的功率譜在區間[0,1/Ts]上近似均勻分布。

時間維持模塊中的迭代式(3)、(4)對輸入序列δtk和輸出序列tk而言構成線性系統。因此,由δtk上的誤差eδtk傳遞到tk的誤差etk也符合線性系統的響應關系,如式(12)所示

etk=etURE,k+etUEE,k

(12)

其中,etURE,k為eδtURE,k經鎖相環傳遞到本地歷元時間估計值tk上的誤差分量。由于eδtURE,k在環路響應的相關時間長度內幾乎為常數,所以,etURE,k≈eδtURE,k。而etUEE,k為eδtUEE,k經鎖相環傳遞到本地歷元時間估計值tk上的誤差分量。對噪聲而言,鎖相環表現為一個帶寬為BL的低通濾波器[19]。即eδtUEE,k的噪聲譜中僅有[0,BL]和[(1/Ts)-BL,1/Ts]兩部分能量能通過環路而到達etUEE,k。因為eδtUEE,k是0均值的,所以etUEE,k也是0均值的。如果BL設計得比較小,以至于相對etURE,k而言,etUEE,k≈0。即etk≈etURE,k≈eδtURE,k。進而,tk的標準差σt即etk的標準差,近似等于eδtURE,k的標準差。根據GNSS接收機誤差傳遞的基本理論,在BL足夠小的條件下,σt的表達式如式(13)所示

σt=TDOP·σURE

(13)

其中,σURE為測距誤差中URE分量的標準差。

2.3 定位誤差值

對誤差傳遞而言,圖2所示算法流程的各模塊均為線性系統。因此,從星上時的觀測誤差到最終的解算誤差之間的傳遞也滿足線性關系。我們可以分別討論URE和UEE對最終定位誤差的影響。

2.3.1 URE導致的定位誤差

觀測量ts中的URE誤差到最終定位誤差的傳播路徑為:

首先,歷史的觀測誤差矢量序列UREm(m

接下來,etk和UREk共同通過方程組(8)傳遞到位置解算結果上,形成由URE導致的誤差分量ePURE,k。

由于各星URE值在幾分鐘的時間內保持不變,可以認為UREm≈UREk。這樣,對于URE誤差分量的傳播而言,方程組(8)與方程組(2)為同解方程。因此,ePURE,k的標準差σP,URE如式(14)所示

σP,URE=PDOP′·σURE

(14)

其中,PDOP′為與G′對應的位置精度因子。

2.3.2 UEE導致的定位誤差

根據2.2節的討論,當時間維持模塊的噪聲帶寬足夠小時,tk中不包含UEE所導致的誤差分量。因此,觀測量ts中的UEE誤差直接通過方程組(8)傳遞到最終定位誤差。

設第k個歷元的定位結果中由UEE導致的誤差分量為ePUEE,k,其標準差σP,UEE如式(15)所示

σP,UEE=PDOP·σUEE

(15)

其中,PDOP為與方程組(8)線性化后的系數矩陣G對應的位置精度因子。

2.3.3 總定位誤差

根據線性系統理論,總定位誤差的瞬時值序列ePk為由URE導致的誤差分量序列ePURE,k和由UEE導致的誤差分量序列ePUEE,k之和如式(16)所示

ePk=ePURE,k+ePUEE,k

(16)

由于URE和UEE兩者的取值是相互獨立的, 所以ePURE ,k和ePUEE ,k兩者的取值也是相互獨立的。本文所提算法的定位誤差ePk的標準差σP可表示為式(17)

=(PDOP′·σURE)2+(PDOP·σUEE)2

(17)

式(17)第1項,由URE導致的定位誤差,與常規導航解算方法完全相同。而第2項,由UEE導致的定位誤差,與PDOP有關。在高軌場景下,UEE取值遠大于URE的值,因此PDOP的大小決定了本算法的性能。

2.3.4 位置精度因子值

PDOP′由方程組(2)線性化后的系數矩陣G′決定。而PDOP由方程組(8)線性化后的系數矩陣G決定。考察這兩個方程組,G′和G存在以下關系:G′=[GN],其中,N為各元素為1的列矢量,其行數與G的行數相等。

PDOP′由G′的權系數矩陣H′對角線的前3項導出。而PDOP由G的權系數矩陣H對角線的3個元素導出。其中

H′=(G′T·G′)-1

(18)

H=(GT·G)-1

(19)

利用實對稱矩陣的特性和柯西不等式可以證明,對于任意的星座構形,一定存在以下關系:PDOP≤PDOP′。即本算法能夠通過消減UEE的影響而改善定位精度。

為了進一步考察本算法對定位精度的改善程度,考慮以下算例。

假設接收機剛好接收到4顆導航衛星的信號, 其中1號衛星距地球最近,其余3顆衛星到1號衛星的夾角相同,均為φ,并且均勻分布在1號星的周圍,如圖3所示。

圖3 PDOP和PDOP’算例模型的空間構型圖Fig.3 Spatial configuration diagram of PDOP and PDOP’ example models

計算φ從1°到179°范圍內對應的PDOP′和PDOP值,如圖4所示。

圖4 PDOP和PDOP’隨φ的變化Fig.4 The value of PDOP and PDOP’ with respect to φ

考察圖4,首先注意到PDOP曲線始終在PDOP′曲線下方,這與之前的分析相符。當φ等于arccos(-1/3)時,PDOP與PDOP′同時取得最小值1.5。

根據引言部分的分析,高軌接收機可見的導航衛星距地心的夾角一般不超過32°,如圖4中“高軌定位場景”所示區間。可見,在這個區間,PDOP僅為PDOP′的1/4～1/10,本算法有顯著的精度優勢。

對應地面接收機而言,可視衛星分布在高度角50°～90°范圍中, 對應圖4中“地面定位場景” 所示區間。可見, 本文提出的算法對地面應用價值不大(并且, 在地面場景中UEE的值也不大)。

3 高軌場景下算法的仿真

通過仿真,進一步驗證了本算法在高軌場景下的定位精度和定位的連續性。

3.1 高軌場景的搭建

(1)用戶軌道數據

從2021年10月26日9時整開始,根據實踐十七號衛星實測導航數據為參考,利用實時更新且準確的BDS導航星座進行仿真。其中高軌GNSS接收機運行軌道的半長軸約為44 379 km。

(2)GNSS星座構形

以同期BDS星座的實際星歷文件作為接收機觀測的GNSS星座。

(3)導航發射信號強度

發射衛星方向圖按照下圖5所示進行設置[20-21]。

圖5 發射衛星天線方向圖設置Fig.5 The Configuration of satellite antenna direction chart

接收機接收信號功率Pr和發射信號功率PT的關系見式(20)

Pr=PT+GT-L+GR

(20)

其中,GT是發射天線增益,GR是接收天線的增益,通常GR=0 dB,L表示導航信號傳輸損耗。仿真系統根據每個歷元上的導航衛星位置和接收機位置分別計算傳輸損耗L和發射天線的方向性增益GT。由PT、GT、L及GR等參數,得到Pr。Pr和接收信號載噪比C/N0的關系為:C/N0=Pr-Nr,Nr是接收機噪聲功率譜密度,取值為-205 dBW/Hz。

仿真中高軌GNSS接收機靈敏度設置為-177 dBW,對應的C/N0的門限值為28。

(4)觀測量誤差模型

第1步:仿真系統按歷元間隔1 s產生無誤差的導航星和接收機之間距離的仿真數據;

第2步:在距離仿真數據上疊加一個OCXO鐘差的實測數據。該數據是以銫鐘為基準,用時間間隔計數器測量獲得的;

第3步:在第2步基礎上添加URE和UEE兩種噪聲。

URE的構造方法:產生一個服從(0,σURE)高斯分布的隨機數分配給各個導航衛星, 在仿真期間各星的 URE 的取值保持不變。σURE取1.0 m。

UEE的構造方法:每個歷元對每顆衛星單獨產生一個服從(0,σDLL)高斯分布的 UEE噪聲。其中的σDLL是根據各星在該歷元的C/N0仿真結果按式(9)計算得到的。

(5)定位誤差和定位連續性統計

接收機模型在每個歷元上選擇C/N0超過閾值的導航衛星的仿真數據按圖2所示流程進行計算。分別統計方程組(2)和方程組(8)的解算結果,作為常規算法和本文提出的算法的結果。

3.2 仿真結果

3.2.1 導航定位服務的可用性

定義導航定位服務的可用性,即高軌場景下定位服務的可用性=高軌場景下定位算法可定位歷元時間/總觀測歷元時間×100%。兩種算法定位服務的可用性情況如表1所示,可用性的改善主要是因為可用衛星數等于3的情況在仿真時段內占據一定比例(圖6中N=3的區間)。

表1 兩種算法定位服務的可用性情況

1) 常規算法。可定位的歷元所占比例為81.17%。2) 高軌時鐘輔助定位算法。新算法可定位的歷元所占比例為90.33%。3) 由1)和2)可見,新算法提高了導航定位服務的可用性。

3.2.2 定位時段精度

利用式(21)評估定位精度值ΔP

(21)

其中,Δkx、Δky、Δkz分別表示第k個歷元接收機三維坐標的誤差值。

兩種算法誤差值ΔP的變化情況如圖6所示。由表2可見:在高軌場景下,本文提出的時鐘輔助定位算法的精度明顯優于常規算法。

圖6 定位歷元時段內定位誤差值變化情況Fig.6 Change of positioning error values during positioning epoch period

可見,傳統算法單歷元解算誤差最大值的數量級為102m;統計后發現:新算法單歷元定位誤差最大值約為15 m左右,新算法單歷元定位誤差值ΔP大幅減小。

兩種算法定位誤差的標準差如表2所示。

表2 定位誤差標準差仿真和理論值

4 結論

本文提出一種高軌衛星導航接收機時鐘輔助定位算法,利用OCXO的穩定性降低了UEE噪聲對定位精度的影響。理論和仿真試驗表明,本文所提的算法特別適用于高軌導航接收機的應用場景,典型的碼偽距定位誤差最大值可從常規方法的150 m左右改善到不超過15 m,并且本算法能夠在可用衛星僅3顆的條件下維持一段時間的定位解算,因此也能在一定程度上改善高軌導航接收機的定位連續性,但對低軌和地面等應用場景,本算法的效果并不顯著。

本文僅從時鐘輔助的角度改善高軌衛星導航接收機的性能,尚不能完全解決高軌衛星導航接收機所面臨的接收信號弱的問題,工程上還需要綜合運用慣性深組合、矢量跟蹤等技術,才能滿足實用化的要求。