羅逸楠，蔡紅濤，高順祖

(武漢大學電子信息學院，武漢 430072)

全球定位系統(tǒng)(global position system，GPS)作為美國國防部研制，能夠適應全天候的無線電導航具有很高的精度[1]、成本費用低、有很好的時效性的特點。而電離層中的總電子含量(total electron content，TEC)的測定對于GPS導航定位會產(chǎn)生一定的影響[2]。除此之外在載波相位測量時，由于整周模糊度和周跳的存在，造成該過程中準確性出現(xiàn)偏差，因此周跳的探測與修復問題在高精度導航定位系統(tǒng)是一個亟待解決的問題[3]。自20世紀80年代以來，已經(jīng)引發(fā)無數(shù)中外學者對此問題的關(guān)注，并提出了TurboEdit法、高次差法、多項式擬合以及小波變換法等解決方法[4]，結(jié)果證明這些傳統(tǒng)周跳探測修復方法在一定程度上取得了較好的效果但都存在其不足之處。TurboEdit法易受到原始數(shù)據(jù)的干擾，無法探測6周以內(nèi)的周跳[4]；多項式擬合法對小周跳的探測很不靈敏[5]；高次差法無法徹底消除噪聲問題[6]；小波變換法則不能探測周跳的真實值[7]。近年來，中國學者提出了新的改進方法。孟令東等[8]提出了將小波變換聯(lián)合偽距相位組合的方法，進一步提高了計算跳變量的數(shù)值穩(wěn)定性；王力等[9]利用歷元間差分觀測值進行參考站探測周跳避免了偽距噪聲的影響；趙松克等[10]利用多項式和多普勒組合的方法改善了探測的誤差性；劉國超等[11]將雙頻載波相位求差法應用其中，有效解決了周跳值解算多解的問題。雖然這些方法都進行了創(chuàng)新和改進，但仍不能完全適應于隨時可能發(fā)生強擾動的電離層環(huán)境中，尤其考慮到周跳現(xiàn)象會對TEC變化產(chǎn)生影響，上述方法仍無法很好解決隨機出現(xiàn)小周跳以及造成TEC不再連續(xù)變化后的探測與修復問題。

針對上述問題，提出了基于相位電離層殘差(phase ionospheric residual，PIR)的組合Kalman濾波算法進行周跳檢測及修復。針對電離層研究TEC變化情況下，采取差分載波相位技術(shù)的(relative total electron content，RTEC)充當實測數(shù)據(jù)的探測量，然后再利用Kalman濾波對檢測到的周跳進行修正。最后通過武漢站和桂林接收機觀測站得到的實測數(shù)據(jù)進行檢驗，進一步驗證本文方法的有效性。

1 基本原理

1.1 周跳及其對TEC的影響

周跳是一種常見的現(xiàn)象，一般多發(fā)生于全球?qū)Ш叫l(wèi)星系統(tǒng)(global navigation satellite system，GNSS)中，由于載波相位測量技術(shù)過程中，衛(wèi)星信號會由于失鎖現(xiàn)象使得整周計數(shù)在某個時間歷元發(fā)生躍遷甚至是中斷消失[12]。因此對于這一現(xiàn)象進行探測及修復研究具有重要意義[13]。

作為GPS當中主要的觀測量，載波相位以及偽距相比而言，前者的測量在精度上有著更好的效果。然而在載波相位測量時，整周模糊度以及周跳現(xiàn)象存在，在一定程度上會對相位測量的準確性造成影響。一般來說，對于沒有出現(xiàn)任何問題的GPS載波相位，其實際觀測值包括整周以及小數(shù)部分二者共同組成。在測量過程中，其中的小數(shù)部分由于載波相位的特性通?？梢暂^為精確地測量出來，但是相位整周數(shù)的變化是由多普勒積分累計計算得到，所以很容易在實際觀測中受到周跳的影響。

對相位的整周數(shù)這一情況，則需通過信號的跟蹤以及整周計數(shù)部分的實際情況來進行考慮。圖1中出現(xiàn)在相位觀測中的一個周跳現(xiàn)象。在500時間歷元時相位由于某種原因出現(xiàn)了一個跳變，造成理論觀測值和實際觀測值出現(xiàn)了一個偏差，而且會影響接下來時間歷元中的觀測現(xiàn)象，這個過程就是出現(xiàn)了一個典型的周跳現(xiàn)象。

圖1 相位觀測值中的周跳

理論上，載波相位觀測值φ可表示為

φ=a+int(φ)+Fr(φ)

(1)

式(1)中：a為初始的整周未知數(shù)；int(φ)為整周計數(shù)；Fr(φ)為不足一周的小數(shù)部分。

在載波相位觀測過程中，對于歷元間的整周計數(shù)如果由于意外發(fā)生了中斷現(xiàn)象，則當調(diào)整好計數(shù)之后會丟失這些數(shù)據(jù)，造成了計數(shù)錯誤，然而對于式(1)中的小數(shù)部分Fr(φ)并不會隨之發(fā)生錯誤，它仍會保持正確變化，周跳現(xiàn)象也就因此而產(chǎn)生。

對于周跳現(xiàn)象發(fā)生在不同情況下有多種因素所致，一般包括如下方面。

(1)由于傳播環(huán)境例如電離層不規(guī)則結(jié)構(gòu)對信號的散射作用或者吸收作用導致的信號強度的深度衰落，以及載波相位的快速起伏。

(2)如果信號的傳播路徑上存在建筑、樹木等障礙物遮，那么將會對信號造成影響進而引起信號的中斷。

(3)再則如果接收機內(nèi)部結(jié)構(gòu)存在有限的動態(tài)范圍，或者鎖相環(huán)帶寬導致的信號失鎖都會引起周跳的產(chǎn)生。

周跳的影響是存在于各方面的，一方面會對觀測數(shù)據(jù)質(zhì)量產(chǎn)生嚴重的干擾；另一方面對于GPS導航、定位以及測量精度等都會存在一定的限制。所以對于周跳的各種問題自始至今是學者們關(guān)注的焦點，特別是對于周跳的探測以及修復尤為關(guān)鍵[14]。同時由于載波相位發(fā)生了突變，計算得到的TEC也會相應地出現(xiàn)隨機跳變，產(chǎn)生嚴重的誤差。圖2為電離層強烈擾動時觀測到的周跳的一個示例，示例中的數(shù)據(jù)選取為2019年3月21日在桂林站對GPS 10號衛(wèi)星觀測所得，圖2中顯示3月21日02：00—06：00由于多次出現(xiàn)了周跳，造成了TEC頻繁的間斷和突跳，而非是正常的連續(xù)變化。

圖2 周跳對衛(wèi)星TEC觀測結(jié)果的影響示例

1.2 差分載波相位測量TEC的方法

當一列初始相位為φ0的波在介質(zhì)中傳播，經(jīng)過時間t傳播了一段距離，此時，該波的相位變?yōu)?/p>

(1)

式(1)中：φ為相位；ω為波的角頻率；c為光速；f為波的頻率；np為相折射指數(shù)；N為電子的電荷量；TECs為傳播路徑上的總電子數(shù)；φ0為波的發(fā)射時初始相位；S為接收機以及衛(wèi)星兩者之間幾何距離。

由于GPS信標L1和L2的載波頻率不同，因此當出現(xiàn)載波頻率不同的兩列波，其相位不能直接進行比較，為了利用雙頻差分相位法獲得TEC，需要進行頻率歸一化，即將不同頻率的載波相位歸一化為共同的基頻fs上，此時可定義差分載波相位為

(2)

式(2)中：m1=f1/fs=154；m2=f2/fs=120，其中，f1為信標L1的載波頻率，f2為信標L2的載波頻率，fs為歸一化后的基頻，fs=11.23 MHz，φ1為信標L1的波相位改變量；φ2為信標L2的波相位改變量。

由式(2)可以計算得到φ1和φ2，同時代入上述定義中可得

(3)

式(3)中：Δφ為差分載波相位；Δφ0為初始相位差，進一步可得

(4)

式(4)中：TEC0為一個常數(shù)，它由初相位以及整周模糊度兩者共同組成。

從式(4)可以看出，利用雙頻差分載波相位測量得到的TEC與實際TEC之間相差一個常數(shù)，這個常數(shù)的值無法精確確定。換句話說，無法經(jīng)過測量得到絕對TEC。因此一般采取差分載波相位技術(shù)用相對總電子含量(relative total electron content，RTEC)表示。在進行這個測量過程中，載波相位的電離層殘差組合充當實測數(shù)據(jù)的探測量，其隨時間歷元的波動即可表示為RTEC的變化。

2 基于Kalman濾波的周跳探測與修復算法

2.1 卡爾曼濾波

L1和L2兩個頻率信號可能會一起出現(xiàn)失鎖以及周跳，同時也可能僅一個頻率出現(xiàn)失鎖和周跳。實際上，相對于L1來說，L2的頻率信號更可能會發(fā)生周跳或者是失鎖現(xiàn)象。但是相對于電離層的研究來說，更關(guān)注TEC的變化，不需要區(qū)分L1和L2信號的周跳值。

因此可以采用基于相位電離層殘差(phase ionospheric residual，PIR)的組合Kalman濾波算法進行周跳檢測以及修復，采用此種方法對存在的周跳RTEC數(shù)據(jù)能夠得到較為理想的處理結(jié)果。

相對于PIR的組合Kalman濾波來說，從本質(zhì)上來講其仍是基于狀態(tài)空間模型，利用遞推的形式對系統(tǒng)的輸入和輸出進行統(tǒng)計的最小方差估計算法[15]。在這個過程不僅需要考慮系統(tǒng)的過程噪聲的統(tǒng)計特性，還不能忽視測量噪聲可能產(chǎn)生的影響，因此在對于多維的非平穩(wěn)的隨機過程中，該算法具有很好的實效性。

首先，假設(shè)系統(tǒng)的狀態(tài)方程為

xk=φk,k+1xk-1+wk-1

(5)

式(5)中：xk為系統(tǒng)在k時刻的n維狀態(tài)向量；φk,k+1為系統(tǒng)的n維狀態(tài)轉(zhuǎn)移矩陣；wk-1為系統(tǒng)在這個過程中所產(chǎn)生的噪聲序列，同時這個序列符合均值為零，滿足協(xié)方差矩陣Qk符合多元正態(tài)分布，即wk～N(0,Qk)，其觀測方程為

yk=Hkxk+vk

(6)

式(6)中：Hk為n維觀測矩陣；vk為均值為零，協(xié)方差矩陣為Rk的服從正態(tài)分布的n維觀測噪聲序列，即vk～N(0,Rk)。

同時，假設(shè)初始狀態(tài)以及每一時刻的噪聲{x0,w1,w2,…,wk,v1,v2,…,vk}均為相互獨立的存在。

將這個過程的觀測量輸入進濾波器中，卡爾曼濾波主要利用觀測過程中噪聲的統(tǒng)計特性[16]，其實際將估計值(即系統(tǒng)的狀態(tài)參數(shù))作為中間濾波器的輸出項。同時，利用時間系統(tǒng)的自動更新聯(lián)合對觀測值的更新算法，對該系統(tǒng)輸入以及輸出端進行比較和聯(lián)系，當系統(tǒng)輸入新的觀測值時，狀態(tài)估計和狀態(tài)估計的誤差隨即進行更新，原則上按照時間更新和測量更新的方程進行分類，其濾波的遞歸方程如下。

首先，按照時間更新方程，可表示為

(7)

(8)

按照測更新后的方程為

(9)

(10)

(11)

式中，Q為系統(tǒng)噪聲；E表示求均值；kk為先驗值和新的測量值之間的相對權(quán)重的一個卡爾曼增益；Yk為測量后更新的觀測值；pk為更新后的誤差協(xié)方差矩陣；R為協(xié)方差矩陣；＾ 表示估計值;～表示第k次迭代的先驗值;通常選取用來衡量xk的估計精度的正定矩陣pk來表示，為xk對應的協(xié)方差矩陣；I為單位矩陣；kk為在先驗值和新的測量值之間的相對權(quán)重的一個卡爾曼增益，且相對估計值在測量值的噪聲比相對較大的情況下，卡爾曼增益使算法對新的觀測值并不是很敏感，反之，如果當估計狀態(tài)相對于觀測值精確度并非精準時，卡爾曼增益則會使得該算法對新觀測的數(shù)據(jù)值更加敏感。如果利用系統(tǒng)方程(或者動態(tài)方程)以及觀測方程估計出所有需要處理的信號，則可以將周跳現(xiàn)象當作是載波相位觀測量的一種噪聲，利用這一算法進行轉(zhuǎn)化，可搭建其周跳探測以及修復模型[17]。

2.2 周跳探測模型

本文模型采取Kalman濾波算法，其原理本質(zhì)上利用最小方差的估計[18]。對于給定的觀測系統(tǒng)，對其輸入或者輸出觀測數(shù)據(jù)進行最優(yōu)的估計。然而在這個過程中，觀測數(shù)據(jù)仍會存在系統(tǒng)噪聲以及干擾，這將對其造成影響，所以這個濾波的過程即為最優(yōu)估計的過程。

在利用差分載波相位測量TEC時，根據(jù)公式計算可知測量得到的TEC與實際TEC的值之間只相差一個常數(shù)。此時，將采取探測量表示為這個殘差組合，相當于觀測RTEC的變化[19]，利用對RTEC的變化的影響檢測周跳的探測和修復效果。

利用做差法對觀測方程進行計算，計算公式公式為

φ=λL1φL1-λL2φL2

(12)

式(12)中：φ為探測量的殘差組分；λL1和λL2分別為信標L1和L2的載波波長；φL1和φL2分別為信標L1和L2的相位。假設(shè)采樣間隔為T，由于PIR以及消去幾何項、鐘差、多徑效應等，只剩下電離層殘差項。此時更高階差分項則將會趨近于0[20]，對于任何一段歷元，對PIR采取用三次多項式模型，同時將φk+1在φk處泰勒展開可得

(13)

將狀態(tài)值記為

Xk=(φk,φ′k,φ″k,φ?k)T

(14)

由式(14)可建立卡爾曼濾波的狀態(tài)模型為

Xk+1=AXk+Guk

(15)

式(15)中：uk為模型誤差；A為狀態(tài)轉(zhuǎn)移矩陣；G為系數(shù)矩陣；根據(jù)φk+1在φk處泰勒展開可分別得到A和G的表達式為

(16)

測量方程為

φk=CXk+zk

(17)

式(17)中：zk為誤差，其方差為σ2；C為系數(shù)矩陣，C=(1,0,0,0)。

記Xk濾波值為X(k|k)，其中k為第k次迭代，濾波誤差方差陣為P(k|k)，Xk+1預報值表示為X(k+1|k)，預報誤差方差矩陣表示為P(k+1|k)，此時濾波的預測方程為

X(k+1|k)=AX(k|k)

(18)

P(k+1|k)=AP(k|k)AT+GQGT

(19)

測量值修正后可表示為

X(k|k)=X(k|k-1)+Mk[φk-CX(k|k-1)]

(20)

式(20)中：Mk為相對權(quán)重的卡爾曼增益，可表示為

Mk=P(k|k-1)CT[R+CP(k|k-1)CT]-1

(21)

(22)

ek=φk-CX(k|k)

(23)

式(23)中：ek為觀測值與預報值之間的殘差值。

2.3 周跳檢測與修復流程

周跳的檢測和修復在實例中，首先設(shè)置濾波器初始值為

X(0|0)=(φ0,φ1-φ0,φ2+φ0-2φ1,φ3+3φ1-3φ2-φ0)T

(24)

(25)

對于檢測周跳主要分為如下3個步驟，關(guān)鍵在于利用Kalman濾波進行修復，流程圖如圖3所示。

Jk為增益矩陣；u為探測閾值，取u=4μ，其中μ為測量值的標準差；Vk預報殘差值

步驟1首先計算Kalman濾波值X(k|k)：①確定濾波啟動的初始條件，此時k=1；②利用包括式(23)之前的方程計算得到預報值、增益矩陣、預報誤差方差陣；③利用計算方程式得到并儲存狀態(tài)濾波值X(k|k)、濾波誤差方差陣、包括觀測到的濾波值φk=CX(k|k)，當k=k+1時，返回②部分進行直至最后一個觀測數(shù)據(jù)出現(xiàn)。

步驟2預測下一時刻狀態(tài)值，同時選取當前的濾波器進行濾波：①固定區(qū)間預測的預測方程為

X(j|k)=Aj-kX(k|k),j>k

(26)

式(26)中：Aj為在j時刻的狀態(tài)轉(zhuǎn)移矩陣；②利用i-1歷元的濾波值X(i-1|i-1)進行預測可得到方程為

X(i+1|i-1)=Al+lX(i-1|i-1),l≥0

(27)

式(27)中：Al為在l時刻的狀態(tài)轉(zhuǎn)移矩陣。

步驟3周跳的檢測與修復，采用預報殘差法：①設(shè)定探測閾值為4μ，與此同時將計算所得到預報殘差與閾值進行比較，比較后可知如果|ek|>4μ，表示在該范圍內(nèi)發(fā)生了周跳現(xiàn)象，得到的周跳值即為預報殘差ek；②進行周跳的消除[22]，將觀測值加上ek；③然后以當前歷元為起點，重新設(shè)置濾波器參數(shù)，重復步驟1。

3 周跳檢測和修復結(jié)果

由上述卡爾曼濾波的周跳探測模型介紹可知，為了驗證Kalman濾波算法的有效性和可靠性，分別進行模擬周跳探測和實測數(shù)據(jù)的周跳修復。

3.1 模擬周跳的探測

首先選取一段“干凈”且沒有跳變的實際觀測數(shù)據(jù)，計算其載波相位電離層殘差組合作為“理論值”，然后利用人為因素在載波相位的隨機歷元中添加隨機大小的跳變，同樣地計算電離層殘差組合，得到含有模擬周跳的觀測值，利用提出的周跳算法對其進行探測和修復，與此同時選擇比較模擬和探測的周跳值的大小，將修復后的觀測值與理論值進行對比，進而分析本文算法的有效性及精度值[23-24]。

選取武漢站的接收機作為第一個站點，該接收機于2018年8月19日觀測到北斗B01衛(wèi)星的數(shù)據(jù)，同時為了盡可能多地去消除多路徑和衛(wèi)星失鎖的影響，將其截至仰角設(shè)置為60°，可以得到十分“干凈”的數(shù)據(jù)。在載波相位L1和L2中隨機加入了模擬周跳，為了突出Kalman濾波算法對小型周跳的探測能力，特意設(shè)定了幾組典型的比較小的周跳，然后使用該算法進行探測和修復，加入周跳的歷元和大小以及探測到的周跳值如表1所示，其中，人為添加的模擬周跳均為整周數(shù)，即1周、2周、3周；而用載波相位的電離層殘差組合充當實測數(shù)據(jù)探測量后，此時探測到的周跳單位為m，可以看出，所有的加入模擬周跳之后的歷元都能成功探測，并且其探測的誤差均小于0.015 m。

從表1可以看出，給出的幾個比較典型的小周跳探測的結(jié)果，實際對應比較大的周跳時，已經(jīng)有更成熟的方法進行更為精確的修復，因此在實際應用中，通常采用的是二級探測的方法，先使用載波相位減去偽距的方式探測并進行修復大于6周的周跳，然后再采用Kalman濾波算法探測較為小的周跳，隨后也進行了一些模擬探測，并統(tǒng)計了探測誤差，轉(zhuǎn)換成TEC單位后的方差為0.005 4 TECu，其中誤差分布如圖4所示。從圖4中的數(shù)據(jù)走勢分布可以看出，總體上基本符合高斯分布。

表1 模擬周跳的探測表

圖4 模擬周跳探測殘差分布

3.2 周跳探測和修復實例

圖5、圖6為兩個實例所使用的數(shù)據(jù)是2019年3月2日在桂林觀測站觀測得到的GPS 5號(圖5)和GPS10號(圖6)衛(wèi)星的數(shù)據(jù)，對此進行周跳探測以及修復?？梢钥闯?，圖5(a)、圖6(a)是由雙頻載波相位實測數(shù)據(jù)計算得到的RTEC數(shù)據(jù)；利用修正后的數(shù)據(jù)對比[圖5(b)、圖6(b)]可以看出，無論是GPS 5號衛(wèi)星還是10號衛(wèi)星，在0～10 000 s歷元這個時間段里，這兩顆衛(wèi)星的RTEC數(shù)據(jù)中都出現(xiàn)了多次短暫的失鎖和明顯的周跳現(xiàn)象，從而使得RTEC變化不連續(xù)，尤其是GPS 5號衛(wèi)星即使在10 000～20 000 s歷元這個時間段里仍然存在很嚴重的數(shù)據(jù)短暫缺失和周跳現(xiàn)象。而圖5(b)中的黑色散點為TEC隨時間的變化，兩條紅線則分別作為根據(jù)計算得到的變化范圍，如果黑色散點落在紅線之內(nèi)表示沒有出現(xiàn)周跳現(xiàn)象。

為了能夠更清楚顯示兩者之間的變化關(guān)系，圖5、圖6中將縱坐標的變化范圍限定在0附近，可以看出，當RTEC變化是平滑的時候，兩條紅線相距很近，而當RTEC開始出現(xiàn)較大的變化時，它們之間的距離相應變寬，使得這種周跳探測及修正方法可以同時適用于電離層平靜條件和有電離層不規(guī)則結(jié)構(gòu)存在的情況。

sTEC為傾斜總電子含量；ΔTEC為RTEC(即相對TEC)的變化率；TECu為TEC的單位，1 TECu≈0.104 m

圖6 模擬周跳的檢測與修正數(shù)據(jù)圖(GPS 10 衛(wèi)星數(shù)據(jù))

在這兩個實例中，這兩顆衛(wèi)星在0～20 000 s歷元之間有兩段時間變化明顯變大，第一段時間0～10 000 s歷元RTEC出現(xiàn)多次周跳，而第二段時間10 000～20 000 s歷元則沒有，相應地在0～10 000 s歷元有很多散點不在這兩條紅線之間，但是10 000～20 000 s歷元則基本都落在了兩條紅線之間。因此利用桂林觀測站觀測得到的實測數(shù)據(jù)很好地證實了本文方法切實可行。

4 結(jié)論

探討了電離層領(lǐng)域的理論研究中周跳現(xiàn)象的存在、產(chǎn)生的原因及可能產(chǎn)生的影響。針對傳統(tǒng)周跳探測修復方法的不足之處，提出基于Kalman濾波的實時TEC周跳探測與修復方法，探測不同的觀測數(shù)據(jù)類型及其組合，利用雙頻差分載波相位測量。同時采取差分載波相位技術(shù)用RTEC代替無法經(jīng)過測量得到絕對TEC，并將修復后的觀測值與理論值進行對比。選取了2018年8月19日武漢站的接收機觀測到的北斗衛(wèi)星的數(shù)據(jù)進行模擬周跳的探測，最后利用2019年3月2日在桂林觀測站觀測得到的實測數(shù)據(jù)進行檢驗，結(jié)果驗證了該方法的有效性。通過分析研究和實驗驗證，得到如下結(jié)論。

(1)周跳的存在會造成TEC出現(xiàn)了頻繁的間斷和突跳，不再連續(xù)變化。對周跳探測以及修復時，即使是小周跳，Kalman濾波算法也有很好的探測效果，實測數(shù)據(jù)驗證其探測的誤差均小于0.015 m。而對于小周跳的探測效果從殘差分布可以看出總體上基本符合高斯分布，更加表明該方法的切實可行。

(2)在進行實際觀測數(shù)據(jù)的對比分析可以看出，由于周跳使得RTEC變化不連續(xù)時，利用Kalman濾波算法可以明顯優(yōu)化和減少此類現(xiàn)象的發(fā)生，無論是在微弱擾動還是劇烈的環(huán)境下進行數(shù)據(jù)處理，都能高效完好的對隨機發(fā)生的周跳現(xiàn)象進行探測識別以及修復工作。

(3)實例試驗是基于桂林站的接收機的數(shù)據(jù)來進行的，并未考慮不同的地區(qū)和更長的時間的情況。同時該算法的精度可能受實際應用場景和觀測數(shù)據(jù)密集程度等多種因素的影響，因此，如果為了達到詳細的分析對比，則需要進行下一步更加完善的實驗設(shè)計。