蔡巍，陳明劍，周舒涵

(中國人民解放軍戰(zhàn)略支援部隊信息工程大學(xué),鄭州 450001)

0 引言

隨著我國北斗衛(wèi)星導(dǎo)航系統(tǒng)(BDS)的發(fā)展，對高精度定位的完整性、連續(xù)性和準(zhǔn)確性的研究變得越來越重要.在BDS 精密定位中，跟蹤BDS 信號的載波相位時，如果在特定的觀察間隔內(nèi)失去鎖相，則所得的觀察結(jié)果與之前和之后的觀察結(jié)果相比，可能會出現(xiàn)不連續(xù)性，這種不連續(xù)性被稱為小周跳，其嚴重影響導(dǎo)航和定位的準(zhǔn)確性和可靠性.周跳檢測和修復(fù)是實現(xiàn)系統(tǒng)高精度定位的有效途徑之一，目前已開發(fā)了許多周跳檢測和修復(fù)的方法.多項式擬合法適用于靜態(tài)數(shù)據(jù)的周跳檢測，但它易受觀測噪聲的隨機性影響，因此對小周跳的檢測并不理想，當(dāng)周跳在多個時期多次發(fā)生時，檢測效果也會變差[1].電離層殘差組合方法可以檢測到小的周跳，但存在檢測結(jié)果多值的問題，無法準(zhǔn)確地檢測和修復(fù)周跳[2].Melbourne-Wübeena(MW)組合方法具有實時功能，但它無法檢測具有相同系數(shù)的組合的周跳[3].為了克服MW 組合方法的缺點，提出了TurboEdit 方法.該方法將MW組合法與無幾何相位組合法相結(jié)合，特別是用于非差和動態(tài)數(shù)據(jù)的周跳檢測.然而，由于偽距的觀測噪聲，無法檢測到小的周跳，上述方法主要適用于單頻和雙頻數(shù)據(jù).隨著載波觀測頻數(shù)的增加，將三頻甚至多頻數(shù)據(jù)結(jié)合的組合系數(shù)進行周跳檢測和修復(fù)已成為一種趨勢.偽距相位組合和無幾何距離組合(GF)是檢測和修復(fù)周跳的典型三頻方法.偽距相位組合方法可以快速檢測周跳，但因含有偽距觀測噪聲的影響，難以檢測到不敏感的小周跳.GF 組合方法可以探測小周跳，但是其組合系數(shù)矩陣的條件數(shù)很大，可能導(dǎo)致病態(tài)方程的發(fā)生.

針對上述方法的不足，本文首先對三頻偽距相位組合采用小波去噪方法降低偽距觀測噪聲的影響[4]，然后將三頻偽距相位組合與三頻GF 組合相結(jié)合，選擇最優(yōu)組合系數(shù)解決偽距相位組合的周跳檢測不敏感和克服病態(tài)方程的問題；最后，通過最小二乘算法和正交三角(QR)分解算法相結(jié)合，固定周跳的最優(yōu)整數(shù)并修復(fù)周跳.

1 偽距相位GF 組合模型及探測周跳原理

BDS 載波相位和偽距觀測方程如下：式中：Pn和φn分別為BDS 第n頻點的偽距和載波相位觀測值，偽距觀測值單位為m，載波相位觀測值單位為周；ρ為衛(wèi)星到接收機之間的幾何距離，單位為m；c為真空中光的速度，單位為m/s；dtr和 dts分別為接收機和衛(wèi)星的鐘差，單位為s；為電離層的延遲系數(shù)，fn為第n個頻點的頻率；I1為頻率f1的對應(yīng)電離層延遲，單位為m;Trop為對流層延遲，單位為m；λn為BDS 第n頻點信號的波長，單位為m；Nn為第n頻點所對應(yīng)的整周模糊度，單位為周；分別為載波相位和偽距觀測值的觀測噪聲.BDS 具有播發(fā)三頻信號功能，因此，選擇B1I、B2I 和B3I 構(gòu)成三頻信號觀測值組合進行實驗，其信號頻率如表1 所示.

表1 BDS 的三頻信號參數(shù)

由式(1)～(2)可推導(dǎo)出三頻偽距和載波相位的表達式：

式中：l、m、n為載波相位觀測值組合的系數(shù)；下標(biāo)a、b、c為三個頻率偽距觀測值組合的系數(shù)；Pabc為偽距組合測量值，單位為m；λlmn為載波相位組合觀測值的波長，單位為m；εφlmn和εPabc分別為組合后的相位和偽距觀測值的觀測噪聲；f1，f2，f3為信號的三種頻率；Nlmn為組合后的整周模糊度，單位為周；φlmn為三頻載波相位組合測量值，單位為周.

采樣率越高，如1 Hz，電離層越平穩(wěn)，相鄰歷元的電離層延遲變化越小.在這種情況下，可以忽略相鄰歷元的電離層延遲影響.用式(3)減去式(4)，可以得到相鄰歷元的差分方程為

式中：Δ代表相鄰歷元間的差分；ΔNlmn表示無幾何偽距組合的周跳估值，單位為周；偽距組合系數(shù)必須滿足a+b+c=1和a2+b2+c2=min.GF 組合可表達為

式中，組合系數(shù)l、m、n必須滿足l+m+n=0和(lλ1)2+(mλ2)2+(nλ3)2=min.對式(11)做歷元間求差，可以得到GF 組合的周跳估值(單位為周)：

假設(shè)BDS 三個頻率的載波相位觀測和偽距觀測相互獨立，根據(jù)文獻[5-6]我們可以設(shè)置載波相位的噪聲標(biāo)準(zhǔn)差δφ1=δφ2=δφ3=δφ=0.01周，偽距觀測值的噪聲標(biāo)準(zhǔn)差δP1=δP2=δP3=δP=0.3m，根據(jù)誤差傳播定律無幾何偽距組合和GF 組合標(biāo)準(zhǔn)差分別為：

當(dāng)組合觀測值的周跳估值絕對值大于其均方根誤差(RMSE)的4 倍時，置信度為99.9%，周跳可以認定為發(fā)生，其RMSE 的4 倍可以視為檢測閥門值[7].當(dāng)發(fā)生周跳時，.

為了更好地檢測和修復(fù)不敏感周跳，選擇了三個線性獨立的組合，包括兩個三頻偽距相位組合和一個三頻GF 組合.根據(jù)文獻[8]，相對于BDS 的B1I、B2I或B3I，無幾何偽距相位組合的載波相位系數(shù)選擇為(0,1,-1)和(1,4,-5)，本文參考文獻[9]中的無幾何相位系數(shù)原則，選擇(-1,1,0)作為GF 組合.

2 小波去噪原理

小波變換在信號局部特征的多尺度分析中具有優(yōu)勢.Daubechies 小波系統(tǒng)的db12 小波變換不僅可以消除噪聲，而且具有較好的正則性和緊支集[10-12]，可以保持信號的階躍或突變點的位置不變.因此，db12 小波變換可用于對GNSS 偽距觀測誤差進行去噪，以提高周跳檢測的性能.

假設(shè)P(t)為偽距信號的離散采樣數(shù)據(jù)，使用db12 小波變換在不同的分辨率下分解，設(shè)C0,k=P(t)，則分解公式為

式中：C為比例系數(shù)；對于偽距信號k為觀測的歷元數(shù)；下標(biāo)j為分解的尺度；D為細節(jié)系數(shù)；n為濾波器長度；h和g分別為高頻和低頻濾波器.考慮到偽距信號的特點，本文選取j=1，即小波分解尺度為一層.

在含有噪聲的原始信號中，無噪聲的有用信號比較平穩(wěn)，一般表現(xiàn)在原始信號的低頻部分，而噪聲信號表現(xiàn)在原始信號的高頻部分.當(dāng)對帶噪聲的原始信號進行小波變換時，小波系數(shù)相對于有用信號的最大范數(shù)將隨著尺度增大，相對于噪聲信號的最大范數(shù)會隨著尺度的增大而變小.因此，有用信號和噪聲信號小波系數(shù)通過小波變換可以清楚地區(qū)分.基于此原則，可以合理設(shè)置高頻系數(shù)的閾值.閾值的選擇和處理策略是小波去噪的關(guān)鍵步驟，Visushrink、SureShrink、Minimax 和BayesShrink閾值選擇和處理策略是常用的方法[13].本文采用Visushrink 閾值確定方法，其表達式為

式中：MMAD為高頻系數(shù)絕對值的中間數(shù)；N為信號長度.閾值消噪的過程中通常采用硬閾值消噪和軟閾值消噪方法.軟閾值函數(shù)具有良好的連續(xù)性和平滑性，但是所有的系數(shù)都被壓縮了，導(dǎo)致信號的一些有效信息丟失.然而，硬閾值函數(shù)不會壓縮信號的階躍點或突變點，因此硬閾值處理策略被廣泛用于去噪.硬閾值函數(shù)可以表示為

式中，Dj,k為小波系數(shù)，經(jīng)過小波變換和去噪后的偽距觀測值表達式為

3 實驗分析

本次實驗采用的是澳大利亞MTIS00AUS0 站點于2022 年1 月30 日采集的數(shù)據(jù)，采樣間隔為1 s，其中BDS 選擇C01 星的B1I、B2I 和B3I 信號作為實驗對象，且信號均不含周跳和粗差.

3.1 無幾何偽距相位小波去噪后性能

根據(jù)第2 節(jié)描述，考慮先使用小波閾值消噪對無幾何偽距相位觀測值進行處理，降低原始偽距觀測值的噪聲水平.圖1(a)、(b)給出了消噪前的無幾何偽距相位組合值，其中圓點表示使用未消噪的無幾何偽距相位組合值得到的超寬巷模糊度N；圖1(c)、(d)給出了消噪后的組合計算得到的超寬巷模糊度N，使用圓點標(biāo)記各個歷元的值.

圖1 不同組合消噪前后得到的超寬巷模糊度N

由圖1 可知：經(jīng)過消噪處理后，混雜在有用信號中的噪聲得到了大幅削弱，其中(0,-1,1)消噪后組合的超寬巷模糊度最大值為14.5188、最小值為14.444 8.而未消噪的超寬巷模糊度最大值為14.524、最小值為14.445 4；(1,4,-5)消噪后組合的超寬巷模糊度最大值為-68.654 2、最小值為-68.918 5，而未消噪的超寬巷模糊度最大值為-68.656 4；最小值為-68.923 7.消噪后波動幅度低于消噪前波動幅度，(0,-1,1)組合消噪前模糊度波動最大為0.078 6，消噪后模糊度波動最大為0.074；(1,4,-5)組合消噪前模糊度波動最大為0.267 3，消噪后模糊度波動最大為0.264 3.

3.2 無幾何偽距相位周跳探測性能

在無周跳的情況下，由圖2(a)、(b)可知，兩個無幾何偽距相位組合中的(0,-1,1)組合周跳探測量的波動均保持在[-0.03,0.03]，(1,4,-5)組合周跳探測量的波動均保持在[-0.15,0.15].由圖2(c)可知，GF 組合(-1,1,0)周跳探測量的波動均保持在[-0.003,0.003].

圖2 不同組合在無周跳下的探測量

由此可知，相較于偽距相位組合，GF 組合的周跳探測能力更高，這是因為GF 組合比無幾何偽距相位組合受到的噪聲影響更小.

原始數(shù)據(jù)不含周跳，為了測試組合的周跳探測能力，本文分別在第50、99、149、199、299、399、499 七個歷元中分別加入(1,0,0)、(0,1,0)、(0,0,1)、(1,1,1)、(5,4,4)、(2,2,2)、(14,14,15)周跳值，這些人為加入的周跳值有小周跳、大周跳、不敏感周跳和特殊周跳，用來驗證本文所選線性組合的周跳探測及修復(fù)能力，圖3 顯示了三種組合觀測模型的探測周跳能力，超出藍色閥值線的紅點為各歷元產(chǎn)生的周跳.

圖3 不同組合在有周跳情況下的探測量

由圖3(a)可知，偽距相位組合中的(0,-1,1)組合可以探測出(0,1,0)、(0,0,1)和(14,14,15)這三種周跳值；由圖3(b)可知，偽距相位組合中的(1,4,-5)組合可以探測出(1,0,0)、(0,1，0)、(0,0,1)、(5,4,4)和(14,14,15)這五種周跳值；由圖3(c)可以看出GF 組合(-1,1,0)可以探測出(1,0,0)、(0,1,0)、(1,1,1)、(5,4,4)、(2,2,2)和(14,14,15)這六種周跳值.綜合起來，將無幾何偽距相位組合與GF 組合結(jié)合起來具有探測小周跳、大周跳、不敏感周跳和特殊周跳的能力，極大地減少了探測盲區(qū).

3.3 周跳探測修復(fù)

當(dāng)載波相位觀測值中的周跳被成功探測出后，需要利用最小二乘算法對周跳進行修復(fù)，通過兩個不同無幾何偽距相位組合的系數(shù)和一個GF 組合的系數(shù)，形成最小二乘中系數(shù)矩陣A，X為三個頻點上待求的周跳浮點解，而L為三個組合在各歷元上的周跳探測量，由式(10)、(12)中的ΔNlmn和ΔNc構(gòu)成.其表達形式如下：

理論上根據(jù)式(19)，可以從周跳的浮點解X直接四舍五入得到周跳的整數(shù)解，由于系數(shù)矩陣中部分元素值較小，周跳探測量L發(fā)生微小的變化會造成方程求解產(chǎn)生較大誤差變化，這就形成了病態(tài)方程組.為了得到更穩(wěn)定的解，本文對系數(shù)矩陣A采取QR 分解，并且在浮點解法確定周跳修復(fù)值初值的基礎(chǔ)上，利用空間搜索和L1 范數(shù)最小原理進一步搜索確定，其表達式為

具體周跳探測量和周跳修復(fù)結(jié)果如表2 所示.

表2 各組合的周跳探測量與GFPMP-GF 組合在去噪前后周跳修復(fù)結(jié)果

由表2 可知，三種組合聯(lián)合進行周跳探測可以覆蓋所有類型的周跳，周跳浮點解和真實的整數(shù)解之間的差值可以用來反映觀測值的噪聲情況，噪聲越大，浮點解與真實的整數(shù)值之間的差值就越大.去噪后，浮點解與真實值之間的差值最大為0.127 9、最小為0.018 4、平均為0.077 0；未去噪時，浮點解與真實值之間的差值最大為0.336 7、最小為0.054 6、平均為0.147 2.由表2 可知，不管去噪還是未去噪，本文提出的算法均能正確地修復(fù)周跳，但通過浮點解和真實值之間的差值比較，采用小波變換去噪后，大大地提高了周跳的探測與修復(fù)性能.

4 結(jié)論

針對BDS 的三頻信號探測與修復(fù)周跳問題，本文提出一種基于小波變換去噪的方法，將無幾何偽距相位與GF 組合模型，并利用真實的北斗數(shù)據(jù)進行實驗驗證后果表明，小波變換去噪的方法可以降低偽距噪聲對無幾何偽距相位組合模型周跳探測和修復(fù)的影響，去噪后的無幾何偽距與GF 相結(jié)合提高了組合模型周跳探測與修復(fù)的性能，可以探測并修復(fù)各類型周跳.