潘毅，張?zhí)祢U，張剛，馬寶澤

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

1 引言

全球?qū)Ш叫l(wèi)星系統(tǒng)（GNSS, global navigation satellite system）能夠提供全天候、連續(xù)且實(shí)時(shí)的高精度地理位置信息及導(dǎo)航、授時(shí)信息、短分組通信等服務(wù)，其中北斗導(dǎo)航定位系統(tǒng)（BDS, Beidou system）是我國獨(dú)立自主研發(fā)建設(shè)的[1-2]。目前，北斗導(dǎo)航定位系統(tǒng)除了采用傳統(tǒng)的二進(jìn)制相移鍵控（BPSK, binary phase shift keying）信號外，提供的L1導(dǎo)航信號等還采用了能夠?qū)崿F(xiàn)頻譜共享及跟蹤性能更優(yōu)越的新信號制式二進(jìn)制偏移載波（BOC,binary offset carrier）信號[3-4]。

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

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

2 高動(dòng)態(tài)BOC信號模型

接收到的高動(dòng)態(tài)BOC信號可以表示為

其中，A是信號幅度；d(t)是導(dǎo)航電文；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是由接收機(jī)和衛(wèi)星之間做徑向運(yùn)動(dòng)產(chǎn)生的多普勒頻偏，kd是由加速度引起的一階多普勒變化率。

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

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

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

其中，Δτ可以表示為

其中，fr表示L1導(dǎo)航信號的射頻頻率。

3 捕獲算法原理

3.1 FRFT算法

時(shí)頻域(t,f)經(jīng)過旋轉(zhuǎn)角α轉(zhuǎn)換到了變換域(u,v)，信號在這個(gè)域的表示則由FRFT給出，則信號x(t)的FRFT定義為

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

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

那么S(n)在p階分?jǐn)?shù)階傅里葉域上的輸出為

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

分?jǐn)?shù)階傅里葉域的能量峰值所在位置可以表示為

其中，是關(guān)于加速度的FRFT旋轉(zhuǎn)角度的估計(jì)值，是關(guān)于頻率誤差的FRFT輸出階數(shù)估計(jì)。多普勒頻偏估計(jì)值多普勒一階變化率估計(jì)值和歸一化FRFT輸出分別為[11]

其中，是旋轉(zhuǎn)角度估計(jì)值那么，關(guān)于多普勒頻偏誤差的FRFT搜索范圍可以表示為

3.2 二階多普勒變化率

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

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

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

其中，ak為實(shí)相位系數(shù)。

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

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

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

將接收到的信號首先經(jīng)過一個(gè)平方環(huán)處理，消除數(shù)據(jù)位和偽碼調(diào)制的影響，只剩下含有高動(dòng)態(tài)項(xiàng)的載波信號，再對其進(jìn)行DPT處理。通過式(19)估計(jì)其階數(shù)，得到的二階DPT幅值譜波形具有如下特征[15]。

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

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

圖2 二階DPT幅值譜波形

4 算法性能和算法流程

4.1 降階處理和FRFT的信噪比

通過瞬時(shí)互相關(guān)的方法對含有高動(dòng)態(tài)項(xiàng)的信號進(jìn)行降階處理，根據(jù)式(1)和式(14)進(jìn)行一次瞬時(shí)互相關(guān)后，信號和噪聲的輸出功率可以分別表示為

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

根據(jù)FRFT的Parseval關(guān)系，其FRFT域的檢測信噪比可以表示為

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

4.2 檢測概率

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

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

非中心參數(shù)λ可以寫成

其中，Tcoh是積分間隔。對于一個(gè)給定的系統(tǒng)虛警概率Pfa，檢測概率為

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

其中，K是懲罰因子，q是搜索單元數(shù)，τd是整個(gè)相關(guān)所需的積分時(shí)間。

4.3 所提算法流程

所提算法步驟如下。

步驟1接收端將接收到的高動(dòng)態(tài)BOC信號下變頻，經(jīng)過模數(shù)轉(zhuǎn)換器（ADC, analog-to-digital converter），以偽碼速率進(jìn)行采樣，得到離散信號r。

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

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

步驟4對本地產(chǎn)生的PN序列進(jìn)行副載波調(diào)制，并根據(jù)階數(shù)確定是否對本地 BOC序列進(jìn)行降階處理，然后與s進(jìn)行相關(guān)處理。

步驟5對相關(guān)結(jié)果進(jìn)行FRFT，根據(jù)式(10)進(jìn)行峰值二維搜索。

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

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

5 仿真實(shí)驗(yàn)

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

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

圖4 Matlab操作界面

5.1 FRFT輸出峰值譜

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

圖5 FRFT輸出峰值譜

5.2 相關(guān)函數(shù)振幅譜

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

圖6 相關(guān)函數(shù)歸一化振幅譜

5.3 檢測概率

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

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

5.4 平均捕獲時(shí)間

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

5.5 計(jì)算復(fù)雜度

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

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

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

對于存在二階多普勒變化率的 BOC信號，所提算法需要進(jìn)行一次瞬時(shí)互相關(guān)運(yùn)算，DPF算法需要進(jìn)行2次瞬時(shí)互相關(guān)運(yùn)算。在實(shí)際應(yīng)用中，當(dāng)時(shí)延長度為，長度為L的BOC序列經(jīng)過一次瞬時(shí)互相關(guān)后，其序列長度變?yōu)椋龠M(jìn)行FRFT，此時(shí)；第二次時(shí)延長度變?yōu)?，進(jìn)行第二次瞬時(shí)互相關(guān)，最后進(jìn)行 PMF-FFT，此時(shí)則所提算法和DPF算法的總計(jì)算量為

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

6 結(jié)束語

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