本發(fā)明涉及無線通信技術(shù)領(lǐng)域，特別涉及一種基于主從式AUKF算法的高動態(tài)衛(wèi)星導(dǎo)航信號載波追蹤方法。

背景技術(shù)：

目前，有關(guān)高動態(tài)接收機載波跟蹤研究可以歸結(jié)為兩個方面：外部其它傳感器輔助載波跟蹤和載波跟蹤算法的改進。

一方面，外部其他傳感器輔助載波跟蹤典型的應(yīng)用是慣導(dǎo)輔助。因為載波環(huán)本身是一種跟蹤較為緊密的環(huán)路，所以，來自外界慣導(dǎo)系統(tǒng)的速度輔助信息必須及時和精確，否則會誤導(dǎo)跟蹤環(huán)路并最終導(dǎo)致信號失鎖。然而，精確的慣性系統(tǒng)通常意味著昂貴的價格，并且外界對其參數(shù)校正的這一過程有可能相當(dāng)復(fù)雜。所以，在沒有慣導(dǎo)輔助時，研究高動態(tài)接收機的載波環(huán)路的跟蹤算法也是當(dāng)今熱點研究課題之一。

另一方面，研究人員提出了多種載波環(huán)路跟蹤算法。有研究者提出了FLL輔助PLL的載波追蹤方法，即通過設(shè)定一個經(jīng)驗閾值讓環(huán)路在兩者之間切換；也有研究者提出了小波變換的載波追蹤方法，即在PLL的鑒相器和環(huán)路濾波器之間加入小波變換，小波變換可以對鑒相器的輸出進行噪聲消除；還有人提出了一種模糊PLL載波跟蹤結(jié)構(gòu)，并利用梯度下降法和遺傳算法優(yōu)化模糊參數(shù)的載波追蹤方法；以及提出了用無跡卡爾曼濾波(UKF)對載波頻率進行估計的跟蹤算法和基于平方根無跡卡爾曼濾波(SR-UKF)的載波跟蹤算法。

載波跟蹤算法中，UKF濾波算法是經(jīng)典卡爾曼濾波(KF)的擴充，它是以無跡變換(UT變換)為基礎(chǔ)，以卡爾曼線性濾波為框架，逼近非線性分布，在一步預(yù)測方程中，使用UT變換來處理均值和協(xié)方差的非線性傳遞的一種非線性濾波方法。相比擴展卡爾曼濾波(EKF)不需要線性化過程，不需要計算雅克比矩陣，具有比EKF更高的估計精度。

但是，UKF算法同KF和EKF—樣，通常先驗假設(shè)過程噪聲和觀測嗓聲為統(tǒng)計特性已知不變的白噪聲，這與高動態(tài)條件下衛(wèi)星導(dǎo)航接收機的工作環(huán)境不符，因此不能保證載波跟蹤系統(tǒng)的收斂性和穩(wěn)定性，即在先驗的噪聲統(tǒng)計特性與實際的噪聲統(tǒng)計特性不相符時，狀態(tài)估計性能將變差甚至發(fā)散。

技術(shù)實現(xiàn)要素：

基于此，針對上述問題，本發(fā)明的目的在于提供一種基于主從式自適應(yīng)無跡卡爾曼濾波算法AUKF的高動態(tài)衛(wèi)星導(dǎo)航信號載波追蹤方法，其能自適應(yīng)調(diào)整過程噪聲方差，從而達到減小模型估計誤差，抑制濾波發(fā)散的目的，使得在高動態(tài)下噪聲統(tǒng)計特性發(fā)生變化時，載波追蹤仍然具有很好的穩(wěn)定性。

本發(fā)明的技術(shù)方案是：一種基于主從式AUKF算法的高動態(tài)衛(wèi)星導(dǎo)航信號載波追蹤方法，包括以下步驟：

a、將接收的衛(wèi)星信號經(jīng)過偽碼剝離后得到信號y，信號y分別與本地壓控振蕩器NCO產(chǎn)生的相位差為90°的兩個復(fù)現(xiàn)載波混頻得到信號I和信號Q；

b、所述信號I和信號Q分別經(jīng)過積分清零電路，得到信號y_I和信號y_Q作為觀測量，并計算信號y_I和信號y_Q的向量形式Y(jié)(k)；

c、所述向量Y(k)通過主從式AUKF濾波算法處理，得到多普勒頻移變化的相關(guān)狀態(tài)向量

d、所述狀態(tài)向量通過環(huán)路濾波電路對頻率估計值進行修正，得到修正值M(k)；

e、用所述M(k)值控制所述本地壓控振蕩器NCO，使其產(chǎn)生的載波頻率和接收信號載波頻率保持一致；

其中，k為離散時間，k＝1，2，3…。

作為對本發(fā)明的進一步改進，所述步驟b中的積分清零電路是以累加清零的方式實現(xiàn)的，每次累加的點數(shù)N為中頻采樣頻率與載波跟蹤環(huán)路的更新頻率之比；當(dāng)完成一次N點累加后，輸出累加結(jié)果，并將累加寄存器清零；累加時鐘為中頻采樣時鐘，清零周期為環(huán)路的更新周期。

作為對本發(fā)明的進一步改進，所述步驟b中，對所述觀測向量Y(k)，有：

其中，L^T＝[1 0 0 0]；X(k)＝[θ(k) ω₀(k) ω₁(k) ω₂(k)]^T，θ(k)為k時刻載波運動參數(shù)，ω₀(k)為k時刻載波多普勒角頻率、ω₁(k)和ω₂(k)為k時刻多普勒角頻率變換率和角頻率變化率的變化率；A為接收信號的幅度，由于不是需要估計的參數(shù)，可設(shè)其為單位值；測量噪聲向量n^T(k)＝[n_I(k) n_Q(k)]為零均值高斯白噪聲，其協(xié)方差矩陣：

其中，n_I(k)表示所述信號I在k時刻的零均值高斯白噪聲，n_Q(k)表示信號Q在k時刻的零均值高斯白噪聲；n_I(l)表示所述信號I在l時刻的零均值高斯白噪聲，n_Q(l)表示信號Q在l時刻的零均值高斯白噪聲；I′為單位矩陣，σ²表示加加速度方差。

作為對本發(fā)明的進一步改進，所述接收信號離散形式下的狀態(tài)向量X(k)是根據(jù)高機動目標跟蹤躍動模型Jerk Model模型確定的，具體過程如下：

b11、定義在Jerk Model模型下，高動態(tài)運動目標加加速度滿足如下指數(shù)衰減規(guī)律：

其中，α為衰減速率，γ為比例常數(shù)；w(t)為t時刻的激勵白噪聲，其方差即策動噪聲強度為為加加速度方差；j(t)為加加速度即加速度在單位時間內(nèi)的變化量，t為連續(xù)時間，t＞0；

b12、根據(jù)所述步驟b11中的衰減規(guī)律，計算載波運動參數(shù)θ(t)分別與載波多普勒角頻率ω₀(t)、多普勒角頻率變換率ω₁(t)和角頻率變化率的變化率ω₂(t)的關(guān)系式：

其中，θ(t)為載波信號中多普勒頻移引起的t時刻相位變化；

b13、定義步驟b12中所述關(guān)系式在連續(xù)時間形式下的向量形式為：

其中,X(t)為連續(xù)形式下t時刻的狀態(tài)向量，表示單位時間內(nèi)X(t)的變化量，A和B為系數(shù)矩陣,并且有：

X(t)＝[θ(t) ω₀(t) ω₁(t) ω₂(t)]^T

B＝[0 0 0 1]^T

b14、根據(jù)步驟b13中所述的向量形式推導(dǎo)離散時間形式下的狀態(tài)向量滿足如下關(guān)系：

X(k+1)＝ΦX(k)+u(k)

其中，Φ為狀態(tài)轉(zhuǎn)移矩陣，u(k)為k時刻的策動噪聲向量，并且有：

u(k)的協(xié)方差矩陣為：

其中，T為采樣時間間隔，Φ(t_k+1-τ)表示(t_k+1-τ)時刻的狀態(tài)轉(zhuǎn)移矩陣，B(τ)表示τ時刻的系數(shù)矩陣，ω(τ)表示τ時刻的激勵白噪聲，q為策動噪聲強度。

作為對本發(fā)明的進一步改進，所述步驟c中的主從式AUKF濾波算法具體包括如下循環(huán)過程：

c1、主UKF算法過程：用k-1時刻的輸出狀態(tài)向量以及步驟c2估計得出的k時刻的噪聲強度來調(diào)整其時間更新式，得到狀態(tài)向量狀態(tài)向量和所述觀測向量Y(k)用于測量更新式中，最后得到k時刻輸出狀態(tài)向量并同時計算新息協(xié)方差ν_k，進入步驟c2；

c2、從UKF算法過程：用k-1時刻的噪聲強度來調(diào)整其時間更新式，得到噪聲強度噪聲強度和所述步驟c1產(chǎn)生的新息協(xié)方差ν_k用于其測量更新式中，最后得到k時刻的噪聲強度再進入步驟c1。

作為對本發(fā)明的進一步改進，所述步驟c1中的時間更新式包括：

χ^*_k|k-1＝f(χ_k-1)

γ_k|k-1＝h(χ_k|k-1)

所述步驟c1中的測量更新式包括：

各參數(shù)存在如下關(guān)系：

λ＝n(ε²-1),

其中，x₀是初始狀態(tài)向量，初始化：P₀為初始狀態(tài)協(xié)方差矩陣，初始化：λ是比例因子；n表示擴維后的系統(tǒng)狀態(tài)向量維數(shù)；ε是散步程度因子，通常取正值；β為驗前分布因子，對于高斯分布，取2最優(yōu)；為一階統(tǒng)計特性時的加權(quán)系數(shù)；為二階統(tǒng)計特性時的加權(quán)系數(shù)；Q^x是先驗過程噪聲協(xié)方差矩陣Q的修正值，其噪聲強度q^x由從UKF估計得到，以逼近Q的噪聲強度；R^x是假設(shè)已知的觀測噪聲協(xié)方差矩陣R。

作為對本發(fā)明的進一步改進，所述步驟c2中的時間更新式包括：

P_k|k-1＝Q^q

所述步驟c2中的測量更新式包括：

各參數(shù)存在如下關(guān)系：

其中，Q^q和R^q分別是過程噪聲協(xié)方差矩陣和測量噪聲協(xié)方差矩陣；為一階統(tǒng)計特性時的加權(quán)系數(shù)；為二階統(tǒng)計特性時的加權(quán)系數(shù)。

作為對本發(fā)明的進一步改進，所述步驟c1中的時間更新式也可為以下簡化式：

所述步驟c1中的測量更新式也可為以下簡化式：

其中，各參數(shù)存在如下關(guān)系：

作為對本發(fā)明的進一步改進，所述步驟d中的環(huán)路濾波電路滿足如下特性：

其中，H(z)是環(huán)路濾波電路的轉(zhuǎn)移函數(shù)，K₀K_d是環(huán)路增益，w_n是自然頻，ξ是阻尼比，T_S是預(yù)檢測積分時間。

本發(fā)明的有益效果是：

(1)觀測向量Y(k)進入主UKF模塊，以主UKF模塊產(chǎn)生的新息協(xié)方差ν_k為依據(jù)，采用從UKF模塊實時估計策動噪聲強度q并自適應(yīng)調(diào)整主UKF的策動噪聲矢量協(xié)方差矩陣Q，使其減小模型估計誤差，抑制濾波發(fā)散；

(2)即便在先驗的噪聲統(tǒng)計特性與實際的噪聲統(tǒng)計特性不相符時，也能通過自適應(yīng)調(diào)整過程噪聲方差，以保持良好的狀態(tài)估計性能，保證載波跟蹤系統(tǒng)的收斂性和穩(wěn)定性；

(3)其能自適應(yīng)調(diào)整過程噪聲方差，從而達到減小模型估計誤差，抑制濾波發(fā)散的目的，使得在高動態(tài)下噪聲統(tǒng)計特性發(fā)生變化時，載波追蹤仍然具有很好的穩(wěn)定性；

(4)通過對主UKF算法過程中各計算式進行合理的簡化，以適應(yīng)在FPGA中的應(yīng)用，極大地提升了運行速率并降低了硬件成本。

附圖說明

圖1是將本發(fā)明實施例所述基于主從式AUKF算法的高動態(tài)衛(wèi)星導(dǎo)航信號載波追蹤方法應(yīng)用于高動態(tài)追蹤環(huán)路的電路原理圖；

圖2是主從式AUKF濾波算法的流程框圖；

圖3是采用傳統(tǒng)UKF濾波算法和本發(fā)明主從式AUKF濾波算法追蹤多普勒頻移對比圖。

具體實施方式

下面結(jié)合附圖對發(fā)明的實施例進行詳細說明。

實施例

參考圖1和圖2，一種基于主從式AUKF算法的高動態(tài)衛(wèi)星導(dǎo)航信號載波追蹤方法，包括以下步驟：

a、將接收的衛(wèi)星信號經(jīng)過偽碼剝離后得到信號y，信號y分別與本地壓控振蕩器NCO產(chǎn)生的相位差為90°的兩個復(fù)現(xiàn)載波混頻得到信號I和信號Q；

b、所述信號I和信號Q分別經(jīng)過積分清零電路，得到信號y_I和信號y_Q作為觀測量，并計算信號y_I和信號y_Q的向量形式Y(jié)(k)；

c、所述向量Y(k)通過主從式AUKF算法處理，得到多普勒頻移變化的相關(guān)狀態(tài)向量

d、所述狀態(tài)向量通過環(huán)路濾波電路對頻率估計值進行修正，得到修正值M(k)；

e、用所述M(k)值控制所述本地壓控振蕩器NCO，使其產(chǎn)生的載波頻率和接收信號載波頻率保持一致；

其中，k為離散時間，k＝1，2，3…。

基于上述方法的高動態(tài)追蹤環(huán)路如圖1，主要包括本地壓控振蕩器NCO、積分清零電路、主從式AUKF算法模塊和環(huán)路濾波電路。

在另一個實施例中，所述步驟b中的積分清零電路是以累加清零的方式實現(xiàn)的，每次累加的點數(shù)N為中頻采樣頻率與載波跟蹤環(huán)路的更新頻率之比；當(dāng)完成一次N點累加后，輸出累加結(jié)果，并將累加寄存器清零；累加時鐘為中頻采樣時鐘，清零周期為環(huán)路的更新周期。

在另一個實施例中，所述步驟b中，對所述觀測向量Y(k)，有：

其中，A為接收信號的幅度，測量噪聲向量n^T(k)＝「n_I(k) n_Q(k)]為零均值高斯白噪聲，其協(xié)方差矩陣：

設(shè)定L^T＝[1 0 0 0],θ(k)＝L^TX(k)，

則測量方程為：

其中，A為接收信號的幅度，由于不是需要估計的參數(shù)，可設(shè)其為單位值；X(k)＝[θ(k) ω₀(k) ω₁(k) ω₂(k)]^T，θ(k)為k時刻載波運動參數(shù)，ω₀(k)為k時刻載波多普勒角頻率、ω₁(k)和ω₂(k)為k時刻多普勒角頻率變換率和角頻率變化率的變化率；n_I(k)表示所述信號I在k時刻的零均值高斯白噪聲，n_Q(k)表示信號Q在k時刻的零均值高斯白噪聲；n_I(l)表示所述信號I在l時刻的零均值高斯白噪聲，n_Q(l)表示信號Q在l時刻的零均值高斯白噪聲；I′為單位矩陣，σ²表示加加速度方差。

在另一個實施例中，所述接收信號離散形式下的狀態(tài)向量X(k)是根據(jù)高機動目標跟蹤躍動模型Jerk Model模型確定的，具體過程如下：

b11、定義在Jerk Model模型下，高動態(tài)運動目標加加速度滿足如下指數(shù)衰減規(guī)律：

其中，α為衰減速率，γ為比例常數(shù)；w(t)為t時刻的激勵白噪聲，其方差即策動噪聲強度為為加加速度方差；j(t)為加加速度即加速度在單位時間內(nèi)的變化量，t為連續(xù)時間，t＞0；

b12、根據(jù)所述步驟b11中的衰減規(guī)律，計算載波運動參數(shù)θ(t)分別與載波多普勒角頻率ω₀(t)、多普勒角頻率變換率ω₁(t)和角頻率變化率的變化率ω₂(t)的關(guān)系式：

其中，θ(t)為載波信號中多普勒頻移引起的t時刻相位變化；

b13、定義步驟b12中所述關(guān)系式在連續(xù)時間形式下的向量形式為：

其中,X(t)為連續(xù)形式下t時刻的狀態(tài)向量，表示單位時間內(nèi)X(t)的變化量，A和B為系數(shù)矩陣,并且有：

X(t)＝[θ(t) ω₀(t) ω₁(t) ω₂(t)]^T

B＝[0 0 0 1]^T

b14、根據(jù)步驟b13中所述的向量形式推導(dǎo)離散時間形式下的狀態(tài)向量滿足如下關(guān)系：

X(k+1)＝ΦX(k)+u(k)

其中，Φ為狀態(tài)轉(zhuǎn)移矩陣，u(k)為k時刻的策動噪聲向量，并且有：

u(k)的協(xié)方差矩陣為：

其中，T為采樣時間間隔，Φ(t_k+1-τ)表示(t_k+1-τ)時刻的狀態(tài)轉(zhuǎn)移矩陣，B(τ)表示τ時刻的系數(shù)矩陣，ω(τ)表示τ時刻的激勵白噪聲，q為策動噪聲強度。

在另一個實施例中，所述步驟c具體包括如下循環(huán)過程：

c1、主UKF算法過程：用k-1時刻的輸出狀態(tài)向量以及步驟c2估計得出的k時刻的噪聲強度來調(diào)整其時間更新式，得到狀態(tài)向量狀態(tài)向量和所述觀測向量Y(k)用于測量更新式中，最后得到k時刻輸出狀態(tài)向量并同時計算新息協(xié)方差ν_k，進入步驟c2；

c2、從UKF算法過程：用k-1時刻的噪聲強度來調(diào)整其時間更新式，得到噪聲強度噪聲強度和所述步驟c1產(chǎn)生的新息協(xié)方差ν_k用于其測量更新式中，最后得到k時刻的噪聲強度再進入步驟c1。

在另一個實施例中，所述步驟c1主UKF算法過程中，Sigma采樣和時間更新的計算公式如下：

χ^*_k|k-1＝f(χ_k-1) (4.2)

γ_k|k-1＝h(χ_k|k-1) (4.6)

所述步驟c1中的測量更新式包括：

各參數(shù)存在如下關(guān)系：

其中，x₀是初始狀態(tài)向量，初始化：P₀為初始狀態(tài)協(xié)方差矩陣，初始化：λ是比例因子；n表示擴維后的系統(tǒng)狀態(tài)向量維數(shù)；ε是散步程度因子，通常取正值；β為驗前分布因子，對于高斯分布，取2最優(yōu)；為一階統(tǒng)計特性時的加權(quán)系數(shù)；為二階統(tǒng)計特性時的加權(quán)系數(shù)；Q^x是先驗過程噪聲協(xié)方差矩陣Q的修正值，其噪聲強度q^x由從UKF估計得到，以逼近Q的噪聲強度；R^x是假設(shè)已知的觀測噪聲協(xié)方差矩陣R。

在另一個實施例中，所述步驟c2從UKF算法過程中，Sigma采樣和時間更新的計算公式如下：

P_k|k-1＝Q^q (4.14)

所述步驟c2中的測量更新式如下：

各參數(shù)存在如下關(guān)系：

其中，Q^q和R^q分別是過程噪聲協(xié)方差矩陣和測量噪聲協(xié)方差矩陣；為一階統(tǒng)計特性時的加權(quán)系數(shù)；為二階統(tǒng)計特性時的加權(quán)系數(shù)。

在另一個實施例中，還提供對步驟c1中的時間更新式和測量更新式進行簡化的方法，以期能實現(xiàn)在FPGA中的應(yīng)用，簡化過程如下：

對χ_k|k-1中的協(xié)方差矩陣進行Cholesky分解，設(shè)定：

聯(lián)合式(4.5)-(4.6)、(7.1)-(7.3)對式(4.7)進行簡化，得到結(jié)果如下：

其中，

a＝W₁+18W₂+2cW₂,

由于，

另外，設(shè)：

將式(7.5)、式(7.6)代入到式(4.8)中，得到觀測值協(xié)方差矩陣如下：

其中，

W₄＝W₆(1-a)²,

W₅＝2W₂,

W₆＝W₃+18W₂,

W₇＝W₁+18W₂,

m＝W₅sin²(c),

n＝2W₂[cos(c)-a]²,

l＝m-n,

關(guān)于狀態(tài)值與觀測值協(xié)方差矩陣的簡化結(jié)果如下：

其中，

綜合上述，整個主UKF濾波算法簡化后可以表示如下，加權(quán)系數(shù)和初始化：

Sigma采樣和時間更新：

測量更新：

由式(7.10)和式(7.11)可以可知，要使能在FPGA中用相應(yīng)模塊實現(xiàn)，主要包括4×4維矩陣的加法、減法，正弦余弦的運算，4×4維矩陣與4×1維矩陣的乘法運算，4×1維矩陣與1×2維矩陣的乘法運算，4×2維矩陣與2×2維矩陣的乘法運算，4×2維矩陣與2×1維矩陣的乘法運算，4×2維矩陣與2×4維矩陣的乘法運算，Cholesky分解等。

在另一個實施例中，所述步驟d中的環(huán)路濾波電路滿足如下特性：

其中，H(z)是環(huán)路濾波電路的轉(zhuǎn)移函數(shù)，K₀K_d是環(huán)路增益，w_n是自然頻，ξ是阻尼比，T_S是預(yù)檢測積分時間。

仿真過程中，采用傳統(tǒng)UKF濾波算法和本發(fā)明主從式AUKF濾波算法追蹤多普勒頻移對比圖如圖3，關(guān)于仿真部分參數(shù)如下。

主UKF中參數(shù)設(shè)置：

載波信號的中頻頻率：l575.42MHz；

采樣時間：T＝lms，載噪比：CNR＝32dB-Hz；

狀態(tài)初始值：x₀＝[0 32995.519 323.356 0]；

初始狀態(tài)協(xié)方差和觀測噪聲協(xié)方差矩陣分別為：

P₀＝diag([0.001 0.001 0.001 0.001])

R＝diag{1×10^-9 1×10^-9}；

UT變換中，取ε＝1.027，κ＝0，β＝2。

從UKF中參數(shù)設(shè)置：

初始噪聲強度：q₀＝5；

初始噪聲強度協(xié)方差：P_q0＝1e-8；

噪聲強度過程協(xié)方差：Q^q＝1e-10；

觀測噪聲協(xié)方差：R^q＝diag{[1e-8 1e-8]}；

UT變換中，取ε＝1.027，κ＝0，β＝2。

以上所述實施例僅表達了發(fā)明的具體實施方式，其描述較為具體和詳細，但并不能因此而理解為對發(fā)明專利范圍的限制。應(yīng)當(dāng)指出的是，對于本領(lǐng)域的普通技術(shù)人員來說，在不脫離發(fā)明構(gòu)思的前提下，還可以做出若干變形和改進，這些都屬于發(fā)明的保護范圍。