GNSS接收機導航濾波器輔助捕獲技術

趙 琳, 羅治斌, 丁繼成, 吳謀炎

(哈爾濱工程大學 自動化學院, 哈爾濱 150001)

全球衛星導航系統(global navigation satellite system, GNSS)的應用已經滲入到了軍用和民用領域[1-2]. 用戶使用GNSS接收機來接收、處理導航衛星發射的電磁信號來確定接收機的位置、速度、時間信息[3-4]. 信號捕獲是接收機使用GNSS導航服務的前提條件，而捕獲裝置的設計需要考慮捕獲性能和資源消耗的平衡. 其中捕獲性能主要指的是靈敏度、捕獲時間，而資源消耗則包括消耗硬件數量和功耗. 在不使用任何先驗信息的前提下，傳統的GNSS接收機信號捕獲算法主要有串行捕獲算法、基于快速傅立葉變換(fast Fourier transform, FFT)的并行頻率捕獲算法、基于FFT的并行碼相位捕獲算法、匹配濾波器(matched filter, MF)捕獲算法、PMF-FFT捕獲算法[5-7]. 串行捕獲算法在載波頻率和碼相位兩個維度進行搜索，硬件消耗較小，但是需要的搜索時間較長；其他的算法屬于并行搜索方法，可以減小搜索時間，但是硬件資源消耗及功耗都會上升.

如果在信號捕獲過程中使用衛星星歷、接收機位置等先驗信息進行輔助，可以提高捕獲性能同時降低資源消耗. 先驗信息可以來自于外部通信網絡，如輔助GNSS(assisted-GNSS)技術；同時也可以來自其他導航系統，如慣性導航系統(inertial navigation system, INS). A-GNSS技術通過外部通信網絡向接收機提供星歷、粗略時間、粗略位置信息來降低搜索區間數量[8-11]. 而接收機位置誤差、動態運動、鐘差、鐘漂增加了搜索的不確定度. 一些研究者研究了INS輔助衛星導航信號捕獲的技術，但是都集中在對衛星多普勒頻率的估計上，在碼相位維度仍然采用串行搜索方式[12-14]. 如果能夠對碼相位估計進行輔助則可以進一步降低搜索區間數量. 由導航濾波器輔助碼相位估計算法可以借鑒矢量跟蹤中對碼環控制量的計算方法，但是矢量跟蹤需要標量跟蹤遷入后才能對碼相位進行計算，無法應用于捕獲階段[15-17]. 目前文獻中并沒有對碼相位直接計算算法的原理、算法的近似處理帶來的影響、誤差源及傳播特性進行深入的分析. 先驗信息同時也可以來自于接收機導航濾波器本身，除了星歷、定位結果、測速結果之外，導航濾波器還可以提供準確的接收機鐘差和接收機鐘漂估計，這些信息都可以提高輔助量的精度，從而減小搜索區間以及捕獲引擎功耗，而導航濾波器輔助捕獲技術缺少相應的理論研究. 因此對GNSS接收機導航濾波器輔助捕獲的研究具有較高的理論價值和實際意義.

本文以北斗B1I信號為研究對象，研究了導航濾波器對信號捕獲的輔助算法、輔助量誤差分析、對捕獲性能的提升. 本文介紹了研究的背景和目前研究存在的不足之處；推導了通過導航濾波器計算偽碼相位、多普勒頻率、NH碼相位的計算方法；分析了導航濾波器輔助碼相位和載波頻率計算的誤差源及傳播特性；簡要介紹了信號捕獲的性能指標；通過理論分析證明了導航濾波器輔助捕獲的優勢.

1 導航濾波器輔助捕獲

1.1 偽碼相位計算

如圖1所示，導航衛星在北斗時(Beidou time, BDT)T1時刻發送的信號在北斗時T2時刻被接收機收到，此時衛星已經由位置1運動到了位置2，其運動矢量定義ΔPsat,而接收機在T2時刻測量的偽距實際上是以T1時刻衛星的位置為基準的. 圖1中其他符號的含義解釋如下：r為衛星與接收機之間的幾何距離，m；τ為信號傳播時間，s；c為光速，m/s；ρ1為T1時刻發送的信號到達接收機所走過的實際路徑長度，m；虛線代表“虛擬路徑”的概念，即假設衛星在位置2發射的信號走過該路徑并在T2時刻到達接收機，虛擬路徑的長度為ρ2. 引入“虛擬路徑”的概念便于對T2時刻導航信號中的偽碼相位進行計算，相當于一條輔助線，而實際上來說信號并不存在這個傳播路徑. 下面將通過“虛擬路徑”的概念來計算T2時刻碼相位的值，計算過程可分為以下4個步驟.

圖1 信號發射時間計算示意圖

步驟1計算虛擬路徑長度. 由于信號傳播時間τ1很短(小于0.2 s)，因此為位置1和位置2之間的距離也很小. 相對于衛星與接收機的距離來說，虛擬路徑和實際路徑幾乎是重合的. 因此對用戶來說，信號在實際路徑和虛擬路徑上的電離層延遲、對流程延遲，以及兩時刻衛星鐘差具有高度的相關性. 由差分定位理論可知，實際路徑和虛擬路徑中電離層延遲、對流層延遲以及兩個時刻的衛星鐘差可以認為是相等的. 虛擬路徑長度可表示為

ρ2=r2+I+T，

(1)

r2=|Psat(T2)-Prec|.

(2)

式中：ρ2為虛擬路徑長度，m；I為電離層延遲，m；T為對流層延遲，m；Psat(T2)為T2時刻衛星在CGCS-2000坐標系內的位置坐標；Prec為T2時刻接收機的位置.

信號在虛擬路徑上的傳播時間可表示為

(3)

式中τ2的值與真實路徑上信號傳播時間τ1近似相等，因此τ2可作為τ1的粗略估計值. 第2節將討論此近似過程所帶來的碼相位計算誤差.

步驟2計算實際傳播路徑長度.T2時刻的測量偽距實際上是以衛星在T1時刻的位置作為基準來定義的. 信號走過的實際路徑長度可定義為

ρ1=r1+I+T.

(4)

式中ρ1為實際路徑長度，m. 衛星在傳播時間τ1(τ1≈τ2)內，由位置1移動到了位置2. 定義Δr為

Δr=r2-r1.

(5)

式中Δr實際上是衛星運動矢量ΔPsat在單位視距矢量上的投影，即

Δr=ΔPsat·l.

(6)

式中l為由接收機指向衛星的單位視距矢量，衛星的運動矢量可以通過兩個時刻衛星位置的差分，得到其計算方法為

ΔPsat=Psat(T2)-C(τ1)Psat(T1)≈

Psat(T2)-C(τ2)Psat(T1)，

(7)

(8)

式中Ωie為地球自轉角速度，rad/s. 對于CGCS-2000坐標系有Ωie=7.292 115 0×10-5rad/s，由于CGCS-2000坐標系屬于旋轉坐標系，方程(7)中C(Δt)矩陣的意義在于將T1時刻衛星位置坐標變換到T2時刻的坐標系中. 此外由于τ1的值不能預知，因此用τ2進行代替. 在通過計算得到Δr后，并且結合式(1)、(4)、(5)可知，實際路徑長度的計算方法為

ρ1=ρ2-Δr.

(9)

信號在實際路徑上的傳播時間τ1可由下式進行計算:

(10)

步驟3計算信號發射時間. 由于接收機在定位之后掌握時間信息，而信號傳播時間已經在步驟2中計算得到，因此在測量時刻處信號發射時間即為接收機時間修正接收機鐘差和衛星鐘差后再減去信號傳播時延，即

ttrans=trec-δtrec+δtsat-τ1.

(11)

式中：ttrans為信號發射時間，s；trec為從接收機時鐘上讀出的T2時刻對應的接收機時間，s；δtrec為接收機鐘差，s；δtsat為衛星鐘差，s. 需要注意的是，trec-δtrec實際上是將接收機時鐘轉換為此時刻的北斗時. 由于信號的發射時間是以衛星鐘作為基準的，因此計算發射時間時還需要將其轉換為衛星鐘上的時間. 繼續補償δtsat這一項將此時刻的北斗時轉換為該時刻衛星鐘時間.

步驟4把信號發射時間轉換為碼相位. 1 ms的發射時間對應一整周的B1I碼(碼相位取值為0～2 046)，而B1I碼的起點與衛星鐘的整數毫秒時間嚴格對齊. 根據BDS B1I信號發射時間與碼相位的關系，可以得到碼相位的計算公式為

φcode(T2)=(ttrans·1 000-[ttrans·1 000])·lcode.

(12)

式中：φcode(T2)為T2觀測時刻碼相位的值，chip；lcode為碼長，chip，對B1I信號來說是2 046；[·]為取整函數;ttrans·1 000實質上將信號發射時間的單位由s轉換為ms.

1.2 載波多普勒計算

由導航濾波器計算多普勒頻率的方法為

(13)

式中：Vsat為T2時刻衛星的速度，m/s；Vrec為T2時刻接收機的速度，m/s；δfrec為接收機時鐘漂移，Hz；δfsat為衛星時鐘漂移，Hz；λ為B1I信號載波的波長，m. 由于衛星上采用了高穩定度的原子鐘，衛星時鐘漂移的量級很小(10-12Hz量級)，因此多普勒頻率計算公式可以簡化為

(14)

1.3 NH碼相位計算

如圖2所示，北斗MEO/IGSO衛星信號中調制有D1電文，其數據位位寬為20 ms. 信號中還調制有NH碼，其周期為20 ms，每個NH碼碼寬為1 ms. NH碼的周期與導航電文數據位持續時間一致. 因此一個導航電文數據比特中包含1組NH碼，也就是20個NH碼符號位[18]. 其造成的影響類似于把20 ms電文寬度變成了1 ms. 頻繁的符號位跳變限制了捕獲過程的相干積分時間，降低了信號捕獲的靈敏度. 在通過導航濾波器估算信號發射時間之后，可以進一步對NH碼相位進行計算. 在獲得NH碼相位后，可以將積分的起始位置選擇為NH碼相位為1的數據塊，這樣便可消除NH碼跳變的影響，從而進行20 ms的相干積分.

圖2 NH碼調制示意圖

NH碼相位與信號發射時間具有嚴格的對應關系. 信號發射時間通常表示為周內秒，取值為[0 s,604 800 s). 信號發射過程中，衛星將發射時間刻在信號上，同時也同步地將NH碼刻在信號上，并且以20 ms為周期不斷循環. 舉例來說，在BDT為0～1 ms間隔中，此時刻入導航信號的NH碼的碼相位為1；而在BDT為19～20 ms間隔中，此時刻入導航信號的NH碼的碼相位為20；在BDT為21～ 22 ms間隔中，由于NH碼的周期性，此時刻入導航信號的NH碼對應相位重新回到1. 因此，對于任意信號發射時刻ttrans，T2時刻信號中的NH碼的碼相位計算方法為

φNH=ceil[rem(ttrans·1 000,20)].

(15)

式中：φNH為T2時刻接收信號中NH碼相位，chip；rem(·)為取余函數；ceil(·)為向上取整函數.

2 輔助算法及輔助量誤差分析

2.1 偽碼相位計算算法誤差分析

在上一小節采用式(7)計算衛星運動矢量ΔPsat時，用傳播時間τ2代替了傳播時間τ1，由此可能會引起ΔPsat的計算誤差，下面對這一項由于近似帶來誤差進行分析. 將式(7)中“≈”兩邊的項進行作差，可得ΔPsat的計算誤差為

δ(ΔPsat)=[C(τ2)-C(τ1)]·Psat(T1)，

(16)

τ2與τ1的時間差為

(17)

式中r1和r2之差主要是由于衛星運動在視距矢量投影所導致的. 北斗GEO/IGSO衛星軌道高度大概為35 786 km，MEO衛星軌道高度為21 528 km. 兩種高度衛星信號傳播到地面上的傳播時間大概為120、72 ms. 另外，GEO衛星速度近似為0，IGSO衛星速度在1～3 km/s之間，MEO速度在2～4 km/s之間. 考慮最差的情況，假定傳播時間τ1=150 ms，衛星與載體相向而行，且相對速度為4 000 m/s，此時r1和r2之差為0.15 s×4 km/s = 600 m. 實際的r1和r2之差應當小于此時的r1和r2之差，有

(18)

將τ1=150 ms，Δτ=2×10-6s代入方程(16)中進行運算，可得

δ(ΔPsat)≈O3×3·Psat(T1)=O3×3.

(19)

因此，即使在最差的情況下，用τ2代替τ1所帶來的ΔPsat計算誤差可以忽略不計.

2.2 偽碼相位計算誤差源及傳播特性

本小節將對導航濾波器計算得到的碼相位及多普勒頻率兩個輔助信息的精度進行分析，尋找輔助量誤差的誤差源和及其傳播特性. 由式(12)可知，碼相位的估計誤差主要取決于信號發射時間的估計誤差. 將式(1)、(7)、(9)代入式(10)，再將式(10)代入式(11)，可以進一步展開T2時刻信號發射時間的計算公式：

ttrans=trec-δtrec+δtsat-

(20)

將公式中的真值用導航濾波器給出的估計值進行替代，可得

(21)

信號發射時間估計誤差可以由估計值減去真值得到，將式(20)和式(21)作差，信號發射時間計算誤差為

(22)

式中：

(23)

(24)

(25)

(26)

(27)

其中：δttrans為發射時間計算誤差，s；δr1為T2時刻由于對接收機位置和衛星位置掌握不準確導致的接收機與衛星之間相對距離估計誤差，m；δI為電離延遲估計誤差，m；δT為對流層延遲估計誤差，m；δ(δtrec)為接收機對鐘差估計的誤差，s，實質上表征了接收機對時間精度的掌握程度；δ(δtsat)為衛星鐘差誤差，s，代表用時鐘模型改正衛星鐘之后殘余的衛星時間誤差，實際表征接收機對衛星時間掌握的精確程度.

下面就式(22)體現的誤差源分別進行討論. 首先對δr1這一項進行分析,如圖3所示，R1為接收機實際位置點，S1為衛星實際位置點，r1為由接收機指向衛星的真實的視距矢量，即

r1=Psat-Prec.

(28)

(29)

(30)

圖3 相對距離誤差與衛星、接收機位置誤差關系

Fig.3 Relation between relative distanceerror and satellite/receiver position error

在直線R1S1上確定A點，使得R1A=R1S2. 則直線AS1的長度即為δr1. 在△R1S2A中，根據余弦定理，有

(31)

δr1=Δr1·l=(δPsat-δPrec)·l.

(32)

將式(32)代入式(22)中，可得

(33)

根據式(33)，可對導航濾波器預測碼的預測誤差的主要誤差源、傳播特性、誤差量級可歸納如下：誤差源主要由衛星軌道誤差δPsat、接收機位置誤差δPrec、電離層延遲估計誤差δI、對流層延遲估計誤差δT，接收機鐘差的估計誤差δ(δtrec)、衛星鐘差的估計誤差δ(δtsat)所構成. 由式(29)可知，δPsat和δPrec經過視距投影后，將會線性地傳播到發射時間計算誤差中去，而其他項則直接以線性的方式傳播到發射時間計算誤差中去.δPrec量級可以參照標準的衛星導航服務的定位精度，其三維均方根大概為7 m，投影到視距上要小于此值，但考慮最差的情況，由此引起的碼相位計算誤差均方根仍可假定為7 m；δ(δtrec)的均方根大概為7 m[19-20]. 衛星位置計算誤差在視距上的均方根大概為2 m；經過模型矯正后δ(δtsat)均方根大概為2 m；經過模型估算后，電離層殘余誤差δI大概在1～5 m之間；經過模型估算后，對流層殘余誤差δT大概在0.1～1 m之間[4]. 綜上所述，碼相位計算誤差的3σ值上限為

5 m+1 m≈36 m.

(34)

2.3 載波多普勒計算誤差源及傳播特性

將式(14)中的真值用導航濾波器的估計值進行替代，可得

(35)

將式(14)與式(35)進行作差，可得到多普勒頻率的估計誤差：

(36)

式中：δVsat為衛星速度估計誤差，m/s；δVrec為用戶速度估計誤差，m/s；δ(δfrec)為接收機鐘漂估計誤差，Hz.

由式(36)可知，δVsat和δVrec經過視距投影后，將線性地傳播到多普勒估計誤差中，而δ(δfrec)將直接線性地傳播到多普勒估計誤差中. 對于BDS B1I提供的導航服務來說，δVrec的三維均方根為0.1 m/s；δVsat的三維均方根應小于0.1 m/s[21]；二者投影到視距矢量后轉化成多普勒估計誤差應小于投影前，但是考慮最差的情況計算多普勒頻率估計誤差3σ值的時候仍然可以不考慮視距投影的影響；接收機時鐘頻率由于受到衛導的校準，其穩定性較高，因此δ(δfrec)量級可假定為(應小于)2 Hz；載波波長λ可以取為B1I信號的波長，大概為0.19 m. 則多普勒頻率估計誤差的3σ值上限為

(37)

由上述分析可知，由導航濾波器計算的碼相位誤差的3σ值大概為36 m(0.24個B1I碼寬)，小于半個B1I偽碼碼寬；而多普勒頻率估算誤差的3σ值大概為6 Hz，小于采用500 Hz為步長時的捕獲頻率精度(250 Hz). 由此可知，采用導航計算輔助捕獲算法，只需要一次相關即可檢測到導航信號是否存在. 需要注意的是，對弱信號采用導航濾波器輔助捕獲需要采用NH碼相位估計和剝離，將相干積分時間擴展為20 ms.

3 捕獲性能指標

3.1 捕獲靈敏度

信號捕獲靈敏度代表滿足一定虛警率和檢測概率的條件下，衛星信號載噪比的最小值. 實際上表征了滿足一定捕獲成功率條件下能夠檢測到的最弱的信號強度. GNSS中頻信號與本地載波和偽碼信號進行混頻、相關運算，得到同相(I)、正交(Q)兩路相關值，其數學模型為

(38)

(39)

(40)

(41)

當信號不存在時，統計量V服從瑞利分布，如下式所示：

(42)

式中fn(v)為瑞利分布的率密度函數.

當信號存在時，統計量V服從萊斯分布，如下式所示：

(43)

式中fs(v)為信號存在時V的概率分布，I0(·)為第一類零階修正貝塞爾函數， 式中a定義為

a=A·R(δτk)·sinc(δfkTcoh).

(44)

若檢測門限選定為Vt，則虛警率Pfa可由下式計算:

(45)

一般情況下，都是先選擇虛警概率(如0.1)，再確定檢測門限. 若首先確定虛警率Pfa，則檢測門限的計算公式為

(46)

在確定檢測門限Vt后，信號的捕獲檢測概率Pd的計算公式為

(47)

在相干積分時間確定的情況下，捕獲概率隨著信號載噪比的上升而上升. 信號捕獲靈敏度定義為:當捕獲概率大于某一數值(如0.95)的條件下，信號載噪比的最小值.

3.2 平均捕獲時間

(48)

4 輔助算法捕獲性能分析

4.1 捕獲靈敏度分析

按照文獻[12]給出的參數設定虛警率、檢測概率、搜索區間和搜索步長：Pfa=0.1，Pd=0.95. 對于無輔助的情況下，其多普勒頻率搜索范圍為[-10 kHz, 10 kHz]，頻率搜索步長為500 Hz，意味著最大的頻率捕獲誤差為250 Hz；碼相位搜索步長為0.5碼片，意味著最大的碼相位捕獲誤差為0.25碼片. 對于有導航濾波器輔助的情況，最大的頻率搜索誤差和碼相位搜索誤差可取為各自估計誤差的3σ值上限：即6 Hz和0.24碼片.

圖4展示了無輔助捕獲情況下和有輔助捕獲情況下對信號檢測概率的分析. 其中輔助算法只采用了導航濾波器輔助碼相位和載波多普勒估計，相干積分時間由于存在NH碼跳變的影響都選擇為1 ms. 從圖中可以看出，當檢測概率Pd=0.95時，無輔助條件下捕獲靈敏度為41.2 dB-Hz，有輔助條件下捕獲靈敏度為40.1 dB-Hz，導航濾波器對捕獲的輔助將捕獲靈敏度提高了1.1 dB.

圖5展示了輔助算法加入對NH碼相位估計的情況下，無輔助和有輔助條件下的捕獲概率分析. 導航濾波器對位同步進行輔助，從而消除了NH碼跳變對相干積分的影響，將捕獲算法的相干積分時間擴展為20 ms；而無輔助的情況下由于NH碼跳變的影響，其相干積分時間仍然被限制為1 ms. 從圖中可以看出，當檢測概率Pd=0.95時，無輔助條件下捕獲靈敏度為41.6 dB-Hz，有輔助條件下捕獲靈敏度為27.3 dB-Hz，導航計算對捕獲的輔助將捕獲靈敏度提高了14.3 dB，相對單純的導航濾波器對多普勒頻率和偽碼相位輔助，NH碼相位輔助可以額外提供12.8 dB的捕獲靈敏度提升.

圖4 有輔助和無輔助條件下的捕獲概率分析

Fig.4 Acquisition probability analysisof aided and unaided acquisition

圖5 增加NH碼相估計條件下捕獲概率分析

Fig.5 Acquisition probability analysis when providing NH code phase estimation

4.2 捕獲時間分析

對于無輔助的二維串行捕獲來說，多普勒搜索區間為±10 kHz，搜索步長為500 Hz，因此頻域搜索區間個數為20 kHz/500 Hz + 1 = 41. 碼相位搜索范圍為[0,2 046)，搜索步長為0.5碼片，因此碼搜索區間數量為2 046/0.5 = 4 092. 因此，二維搜索總共的搜索區間數量為4 092×41 = 167 772. 文中無輔助串行捕獲算法相關參數為:Pd=0.95，Pfa=0.1，K=1，Tcoh=1 ms，M=1，Lc=1. 將這些參數代入式(48)，得二維搜索的平均搜索時間為93.64 s.

對于導航計算輔助捕獲算法，總共需要搜索的區間數量為1，其他參數與無輔助的二維捕獲相同，通過式(48)計算可得需要的平均搜索時間，當Tcoh=1 ms時候平均捕獲時間為1.1 ms；當Tcoh=20 ms時，平均捕獲時間為21.1 ms. 可見在導航濾波器的輔助下，捕獲所需的時間大大降低了，擴展相關積分時間并沒有明顯的提高捕獲時間.

5 結 論

1)在接收機已經完成導航計算，并且掌握導航星歷的前提下，可以采用導航濾波器對信號捕獲進行輔助，從而提高捕獲靈敏度，降低捕獲時間，同時可以降低無輔助條件下捕獲引擎頻繁運行帶來的功耗，更好地平衡捕獲裝置性能與資源消耗的關系.

2)仿真表明，采用導航計算輔助載波多普勒和偽碼相位搜索，可以提高1.1 dB的捕獲靈敏度，平均捕獲時間降低為1 ms；在輔助量包含位NH碼相位的前提下，通過擴展相干積分時間可以進一步擴展相干積分時間為20 ms,從而提供額外的12.8 dB的捕獲靈敏度，平均捕獲時間為20 ms左右.

3)矢量跟蹤環路通過導航計算與基帶跟蹤的融合來實現星間輔助與跟蹤靈敏度提升，所提出的導航濾波器輔助捕獲算法進一步擴展了矢量跟蹤的概念：將導航濾波器先驗信息用來輔助基帶捕獲和位同步，在更廣泛的意義上實現接收機導航計算與基帶處理的融合與相互增強，因此導航濾波器輔助捕獲算法具有一定理論價值和實際價值.