高動態BOC信號捕獲算法

潘毅，張天騏，張剛，馬寶澤

（重慶郵電大學通信與信息工程學院，重慶 400065）

1 引言

全球導航衛星系統（GNSS, global navigation satellite system）能夠提供全天候、連續且實時的高精度地理位置信息及導航、授時信息、短分組通信等服務，其中北斗導航定位系統（BDS, Beidou system）是我國獨立自主研發建設的[1-2]。目前，北斗導航定位系統除了采用傳統的二進制相移鍵控（BPSK, binary phase shift keying）信號外，提供的L1導航信號等還采用了能夠實現頻譜共享及跟蹤性能更優越的新信號制式二進制偏移載波（BOC,binary offset carrier）信號[3-4]。

實際的通信環境日益復雜，飛行器速度快、加速度大，甚至加加速度大，接收機往往是在高動態環境下工作，使接收到的信號具有大多普勒頻偏和一階、二階多普勒變化率。這對BOC信號的捕獲提出了較大的挑戰，而大部分已經提出的捕獲算法[5-7]只能在載波頻率完全消除的前提下使用，不能對加速度和加加速度分別引起的一階、二階多普勒變化率分量進行有效的補償。為了解決高動態下的捕獲問題，文獻[8]采用部分匹配濾波器（PMF, partial matched filter）和快速傅里葉變換（FFT, fast Fourier transform）相結合的算法，提高了多普勒頻偏估計的范圍和分辨率，但是無法對多普勒變化率進行補償。文獻[9]通過離散chirp傅里葉變換對多普勒頻偏和一階多普勒變化率進行估計，但在高速徑向運動中效果不盡如人意。文獻[10]采用最大似然算法來估計碼相位、多普勒頻偏和載波相位誤差，但同時也引入了巨大的計算量，難以在要求實時性、結構簡易的接收機中使用。當偽碼對齊時，擴頻信號的相關結果可以看作線性調頻（LFM, linear frequency modulation）信號，文獻[11-12]引入了分數階傅里葉變換（FRFT, fractional Fourier transform）算法來解決高動態GNSS信號的捕獲問題。

目前，幾乎所有的高動態捕獲算法都是針對大多普勒頻偏或存在一階多普勒變化率的情況，因此本文采用離散多項式相位變換（DPT, discrete polynomial-phase transform）算法來處理存在二階多普勒變化率的BOC信號，并通過先對BOC信號定階再根據信號動態階數分別進行處理，最后通過FRFT算法進行捕獲，最終形成了一套完整的高動態BOC信號捕獲結構。

2 高動態BOC信號模型

接收到的高動態BOC信號可以表示為

其中，A是信號幅度；d(t)是導航電文；s(t)=c(t)sc(t)是基帶信號，c(t)是偽碼（PN, pseudo-noise）信號，sc(t)=sgn[sin(2πfsct)]或 sc(t)= sgn[cos(2πfsct)]是副載波信號，fsc是副載波頻率，φ0是初始相位，n(t)是均值為 0、方差為σ2的加性高斯白噪聲。載波信號φ(t)可以表示為

其中，fi是中頻頻率，fd是由接收機和衛星之間做徑向運動產生的多普勒頻偏，kd是由加速度引起的一階多普勒變化率。

經過相關和積分處理后，載波信號φ(t)的輸出可以表示為

其中，T表示相干積分時間，Δφ0表示載波相位誤差，R(·)表示BOC信號的自相關函數，Δτ表示碼相位時延誤差。假設多普勒變化率在估計過程中不變，則可以得到一個平均角速度為

那么，式(3)可以改寫為[13]

其中，Δτ可以表示為

其中，fr表示L1導航信號的射頻頻率。

3 捕獲算法原理

3.1 FRFT算法

時頻域(t,f)經過旋轉角α轉換到了變換域(u,v)，信號在這個域的表示則由FRFT給出，則信號x(t)的FRFT定義為

其中，k為整數，為旋轉角度，p為FRFT的階數，Fp[·]為FRFT的算子符號。最初的時頻域(t,f)經過旋轉角α轉換到了變換域(u,v)，能夠使原本在(t,f)平面發散的頻譜在(u,v)平面得到聚集，如圖 1所示，從而能夠估計出多普勒頻偏fd和多普勒一階變化率kd。

圖1 FRFT：(t,f)平面變換到(u,v)平面

那么S(n)在p階分數階傅里葉域上的輸出為

其中，S(n)是一個帶噪的線性調頻信號，ω=UΔt，Δt是變換域的采樣間隔。離散信號S(n)關于α的N點離散FRFT可以表示為

分數階傅里葉域的能量峰值所在位置可以表示為

其中，是關于加速度的FRFT旋轉角度的估計值，是關于頻率誤差的FRFT輸出階數估計。多普勒頻偏估計值多普勒一階變化率估計值和歸一化FRFT輸出分別為[11]

其中，是旋轉角度估計值那么，關于多普勒頻偏誤差的FRFT搜索范圍可以表示為

3.2 二階多普勒變化率

r(n)的M階瞬態矩定義為

其中，其中，τ是時延長度。那么，r(n)的M階DPT可以表示為[14]

其中，N1為采樣長度，Ts為采樣間隔，ω為角頻率。M階多項式相位信號可以表示為

其中，ak為實相位系數。

將式(18)代入式(17)，可以得到

當載波存在二階多普勒變化率時，采樣后的BOC信號r(n)可以表示為

其中，kf表示加加速度引起的二階多普勒變化率。通常情況下，接收機和衛星徑向加速度為-100g～100g，所引起的最大一階多普勒變化率為±5.15 kHz/s，而加加速度的范圍是-40g/s～40g/s，其中，g為重力加速度。即最大二階多普勒變化率大約為±2.06 kHz/s2。當kf≠0時，FRFT算法不適用于具有加加速度的信號，所以本文采用瞬時互相關的方法進行降階。由瞬態矩性質可得，P2[s(n),τ]可認為是原 BOC基帶信號s(n)經過瞬時互相關運算變成新BOC基帶信號s′(n)，因此可以得到

將接收到的信號首先經過一個平方環處理，消除數據位和偽碼調制的影響，只剩下含有高動態項的載波信號，再對其進行DPT處理。通過式(19)估計其階數，得到的二階DPT幅值譜波形具有如下特征[15]。

1) 如果二階DPT幅值譜是具有一定帶寬的凌亂譜線，如圖2(a)所示，說明載波信號是三階多項式相位信號并含有二階多普勒變化率，則需要對信號先進行降階處理再采用FRFT估計其動態參數。

2) 如果二階DPT幅值譜是非0頻的單線譜或者只有0頻的譜線，如圖2(b)或圖2(c)所示，說明載波信號分別是二階或一階多項式相位信號，并不含有二階多普勒變化率，可以直接采用FRFT算法對接收到的信號參數進行估計（如果信號只含有多普勒頻偏時，則FRFT算法的階數為1，變為傳統的傅里葉變換）。

圖2 二階DPT幅值譜波形

4 算法性能和算法流程

4.1 降階處理和FRFT的信噪比

通過瞬時互相關的方法對含有高動態項的信號進行降階處理，根據式(1)和式(14)進行一次瞬時互相關后，信號和噪聲的輸出功率可以分別表示為

其中，L為采樣后信號的長度。結合式(24)和式(25)可得輸出信噪比為

根據FRFT的Parseval關系，其FRFT域的檢測信噪比可以表示為

其中，var表示方差，S(α,U)表示信號在FRFT域的幅值，W(α,U)表示噪聲在 FRFT的幅值。FRFT算法能夠使降階處理后的信號在最佳旋轉角度的變換域中能量得到聚集，即var[S(α,U)]取得最大值。但對于高斯白噪聲在任意階次的分數階傅里葉域內能量分布是均勻的，即var[W(α,U)]是不變的，所以在估計動態參數的同時抑制了噪聲并提高了信號的檢測信噪比，增強了系統的抗噪性能，更易于實現高動態信號的快速捕獲。另一方面，對于存在二階多普勒變化率的BOC信號，所提算法只需要進行一次瞬時互相關就可以達到降階的目的，而傳統的采用DPT和PMF-FFT相結合的捕獲算法（以下簡稱為DPF算法）[17]需要進行2次瞬時互相關才可以進行捕獲。可見，在采用相同的延時長度的情況下，所提算法的輸出信噪比損失更小，降階計算量更小。

4.2 檢測概率

由于載波信號中的多普勒變化率和多普勒頻偏能通過FRFT精確地被估計，信號的高動態項不會影響FRFT輸出峰值的檢測，因此基于FRFT捕獲算法的檢測概率和基于 FFT捕獲算法的檢測概率是相同的。接收信號和本地信號相關后的輸出服從 Rice分布，其虛警概率Pfa和檢測概率Pd分別表示為

其中，Pn(x)是只含有噪聲時的概率密度函數，Ps(x)是含有有用信號和噪聲時的概率密度函數，I0是 0階修正貝塞爾函數，V是門限值，可以表示為

非中心參數λ可以寫成

其中，Tcoh是積分間隔。對于一個給定的系統虛警概率Pfa，檢測概率為

其中，Q(x,y)是MarcumQ函數。平均捕獲時間也是捕獲算法一個重要的指標，可以表示為[18]

其中，K是懲罰因子，q是搜索單元數，τd是整個相關所需的積分時間。

4.3 所提算法流程

所提算法步驟如下。

步驟1接收端將接收到的高動態BOC信號下變頻，經過模數轉換器（ADC, analog-to-digital converter），以偽碼速率進行采樣，得到離散信號r。

步驟 2采樣后的信號通過本地載波數字控制振蕩器（NCO, numerically controlled oscillator）進行載波剝離。

步驟3對仍存在殘余頻偏的信號進行定階運算，根據階數確定是否進行降階處理，得到的信號為s。

步驟4對本地產生的PN序列進行副載波調制，并根據階數確定是否對本地 BOC序列進行降階處理，然后與s進行相關處理。

步驟5對相關結果進行FRFT，根據式(10)進行峰值二維搜索。

步驟6將相關結果的最大值與門限V進行比較（V值由式(30)得到）。如果最大值超過預設門限值V，則說明本地PN序列與接收BOC信號的偽碼相位對齊，進入跟蹤階段（本文主要解決捕獲問題，跟蹤過程不予以討論）；如果沒有超過預設門限值V，則通過邏輯控制本地載波 NCO，根據式(11)得到的動態參數對多普勒頻偏進行補償，同時滑動本地PN相位。

步驟7重復步驟3～步驟6，直到偽碼相位對齊。如此就形成了一整套捕獲方案，如圖3所示。

5 仿真實驗

為了檢測所提捕獲算法的性能，本文采用Matlab對 sinBOC(1,1)信號進行仿真實驗，其中信號的采樣頻率為16.368 MHz，碼速率為1.023 MHz，碼片長度為1 024。FRFT運算時間為51.2 ms。仿真參數設置為多普勒頻率fd=50 kHz，一階多普勒變化率為 5 kHz/s，二階多普勒變化率為 2 kHz/s2，SNR=-5 dB。實驗中采用DPF算法進行對比仿真，其中PMF時間為16 ms，碼片長度取3 069，分為3段，每段長度為1 023。圖4為Matlab實驗操作界面。

圖3 本文所提算法捕獲流程

圖4 Matlab操作界面

5.1 FRFT輸出峰值譜

由于高斯函數的Wigner-Ville分布仍然是高斯函數，高斯白噪聲在任何階數FRFT下仍然是高斯形狀。而載波信號能夠在相應階數的FRFT下實現能量聚集，比較易于信號的檢測與參數估計。FRFT輸出峰值譜如圖5所示。由圖5可以明顯看出，輸出譜在p=1.051和u=444處出現了一個峰值。根據式(11)和式(23)可以計算出，多普勒頻偏fd=51.32 kHz，多普勒一階變化率kd=4 930.00 Hz/s，多普勒二階變化率kf=2 008.98 Hz/s2，這與之前的設定值誤差保持在2%左右，因此所提算法能夠準確估計高動態BOC信號的動態參數。

圖5 FRFT輸出峰值譜

5.2 相關函數振幅譜

圖6是對存在二階多普勒變化率的BOC信號分別采用本文所提算法、DPF算法和PMF-FFT算法進行捕獲得到的相關函數歸一化振幅值。圖 6(a)中，所提算法在碼相位為596位置的主峰幅值明顯高于其他位置；圖6(b)中，DPF算法的振幅譜能夠在相同的碼相位處得到一個最大幅值，但有多個位置的幅值與最大幅值非常接近，這就造成接收機無法正確鎖定主峰位置；圖6(c)中，PMF-FFT算法輸出最大幅值所對應的碼相位位置錯誤。通過圖6對比分析，可以證明本文所提算法的有效性，也使判決器更容易檢測到主峰位置，且檢測概率更高。

圖6 相關函數歸一化振幅譜

5.3 檢測概率

對不同動態參數下的 BOC信號采用本文所提算法、DPF算法和PMF-FFT算法得到的檢測概率如圖7所示。從圖7可以看到，對于存在二階多普勒變化率的信號，所提算法比 DPF算法提高了3.5 dB；對于只有一階多普勒變化率的信號，所提算法比 DPF算法提高了 4.5 dB。而傳統 PMF-FFT算法在大動態參數條件下已經無法實現信號的捕獲。

圖7 本文所提算法、DPF算法和PMF-FFT算法的檢測概率

5.4 平均捕獲時間

捕獲時間的參數設置為搜索單元數q=2 048，虛警概率Pfa=0.01，積分時間τd=1 ms，懲罰因子K=1。從圖8可以看出，對于存在二階多普勒變化率的信號，當SNR=-5.5 dB時，本文所提算法的平均捕獲時間降到最低，大約是DPF算法平均捕獲時間的；對于只有一階多普勒變化率的信號，當SNR=-13 dB時，所提算法的平均捕獲時間降到最低，大約是DPF算法平均捕獲時間的。在低信噪比下，所提算法對多普勒頻偏、多普勒變化率等高動態參數的估計更加準確，殘余頻偏較小甚至為0。這使捕獲檢測函數引入的誤差較小，更有利于碼相位的搜索。但是隨著信噪比的升高，DPF算法的平均捕獲時間要低于所提算法，這是由于此時多普勒頻偏估計更準確，平均捕獲時間主要取決于算法的計算時間，DPF算法估計時間為16 ms，FRFT估計時間為51.2 ms[12]，即DPF算法的計算復雜度要低于FRFT算法。而接收到的BOC信號的信噪比往往都是較低的，所以對于算法來說在低信噪比條件下的性能表現更加重要。

5.5 計算復雜度

FRFT算法和 PMF-FFT算法的計算量可以表示為[13]

圖8 本文所提算法和DPF算法的平均捕獲時間

其中，pu是 FRFT階數搜索的上限值，pl是 FRFT階數搜索的下限值；；τ1和τ2分別是FRFT和PMF-FFT算法全部估計所需要的時間；T1和T2分別是FRFT和PMF-FFT算法單個輸出的時間是FRFT和FFT算法輸出的比值。而進行一次瞬時互相關運算需要的計算量為

對于存在二階多普勒變化率的 BOC信號，所提算法需要進行一次瞬時互相關運算，DPF算法需要進行2次瞬時互相關運算。在實際應用中，當時延長度為，長度為L的BOC序列經過一次瞬時互相關后，其序列長度變為，再進行FRFT，此時；第二次時延長度變為，進行第二次瞬時互相關，最后進行 PMF-FFT，此時則所提算法和DPF算法的總計算量為

可見，當K12=1時，即2種算法輸出信號長度相同時，所提算法的計算量比DPF算法要高，這是由于DPF算法的硬件需求更高，在增加匹配濾波器個數的情況下使計算量降低。

6 結束語

對存在二階多普勒變化率的高動態BOC信號，本文引入DPT算法進行先定階后降階的處理，然后基于 FRFT算法估計其多普勒頻偏和碼相位。所提算法的檢測概率、平均捕獲時間和計算復雜度都進行了分析，并在仿真中與DPF、PMF-FFT算法進行了性能對比。對于存在二階多普勒變化率的高動態BOC信號，所提捕獲算法的檢測概率比DPF算法提高了3.5 dB，并在低信噪比情況平均捕獲時間較少。本文提出的一整套捕獲方案不但解決了傳統捕獲算法無法估算二階多普勒變化率的問題，而且提高了系統的抗噪性能，非常適用于高動態BOC信號。