蔡成林，沈文波，曾武陵，于洪剛，謝小平

1. 湘潭大學(xué)自動化與電子信息學(xué)院，湖南 湘潭 411105； 2. 桂林電子科技大學(xué)信息與通信學(xué)院，廣西 桂林 541004； 3. 廣西電網(wǎng)有限公司柳州供電局，廣西 柳州 545005； 4. 廣州大學(xué)經(jīng)濟與統(tǒng)計學(xué)院，廣東 廣州 510006

周跳探測與修復(fù)是載波相位定位中關(guān)鍵的一步，能否準確探測與修復(fù)，直接影響到定位精度。周跳可以是從1到百萬的整數(shù)[1]，大周跳很容易被探測出來，探測最小的1周周跳是精密定位的要求。目前周跳探測的方法有很多，對單頻數(shù)據(jù)而言，簡單的方法有高次差法[2]和多項式擬合法[3]，這兩種方法都能容易地探測出大周跳，但難以探測出小周跳。對于雙頻數(shù)據(jù)而言，常用的是TurboEdit[4]方法，它是聯(lián)合使用MW法和電離層殘差法[5]來探測周跳。對于三頻數(shù)據(jù)的周跳探測方法，多數(shù)學(xué)者采用不同的無電離層組合和MW組合進行求解周跳[6-10]，主要是研究如何選取更優(yōu)的組合系數(shù)。

雖然三頻組合的數(shù)據(jù)可以提供更多更優(yōu)的組合方式[7-10]，但是雙頻數(shù)據(jù)組合方式有限，在低采樣率條件下，要達到高精度探測較為困難。為了更好地推進雙頻探測周跳的方法，文獻[11]提出了在MW基礎(chǔ)上向前和向后移動窗口濾波的算法(FBMWA)和電離層殘差的二次時間差法(STPIR)。其研究表明FBMWA法比文獻[4]提出的用遞歸的均值濾波器對MW進行濾波效果更好，但需要用到后面歷元的數(shù)據(jù)，不利于實時探測周跳。STPIR法與電離層殘差法相比，在電離層變化緩慢時，效果相當(dāng)；在電離層活躍時，電離層殘差法效果較差，而STPIR法依然較好。

另外，文獻[12]首先提出了多普勒積分法探測周跳，該方法簡單易行，但是多普勒的精度受采樣率的約束，即使在1 s采樣率時也未必能保證很好地探測出1周周跳。本文提出了一種重構(gòu)多普勒的方法，重構(gòu)出的多普勒受采樣率的影響小。在低采樣率時，重構(gòu)相應(yīng)頻率的多普勒替代其觀測值，代入多普勒積分法進行周跳探測，效果更優(yōu)。受寬巷組合的啟發(fā)，本文提出了一種組合雙頻載波，并重構(gòu)其多普勒用于積分，進行周跳探測的方法。它只采用相位進行組合，而MW法組合時引入了偽距，偽距噪聲較大，所以它比MW法探測精度更優(yōu)，噪聲更低。且它與STPIR法聯(lián)合求周跳，可以探測出最小的1周周跳以及連續(xù)周跳。

1 重構(gòu)多普勒積分聯(lián)合STPIR

1.1 重構(gòu)多普勒的方法

本文利用衛(wèi)星與接收機間的多普勒形成的原理和模型[13]，如圖1所示，提出了一種重構(gòu)多普勒的方法，僅需要廣播星歷，適合實時應(yīng)用。下面闡述重構(gòu)多普勒的算法模型：由于廣播星歷計算得到的衛(wèi)星位置Ps以及衛(wèi)星速度vs三維向量都是基于地心地固(ECEF)坐標系的，即該坐標系隨地球自轉(zhuǎn)而轉(zhuǎn)動，因此不需要考慮地球自轉(zhuǎn)對多普勒的影響。

圖1 衛(wèi)星與接收機的多普勒模型Fig.1 Doppler model of satellite and receiver

如圖1所示，向量vs表示接收機的速度，衛(wèi)星在接收機處的觀測向量ρ的方向由接收機指向衛(wèi)星。設(shè)vs=(xvs,yvs,zvs)，Ps=(xs,ys,zs)，偽距定位方式得到的接收機位置Pu=-(xu,yu,zu)。故ρ=(xρ,yρ,zρ)=(xs-xu,ys-yu,zs-zu)。那么多普勒重構(gòu)值可由式(1)計算

(1)

式中，l是單位觀測向量，l=ρ/|ρ|；λ是對應(yīng)衛(wèi)星頻率的波長；v是接收機速度向量；“·”表示向量點乘。

1.2 重構(gòu)多普勒積分法

分析式(1)可知重構(gòu)多普勒誤差主要來源于中分子計算的部分，若組合出來的波長越長，即式(1)中分母變大，重構(gòu)誤差就會減小。下面以北斗為例，北斗的頻率B1=1 561.098 MHz和B2=1 207.14 MHz，對應(yīng)的相位觀測值φ1和φ2可以表示成如下[14]

φ1=N1+φ1+c(dtr-dts)/λ1+ΔT+ΔI1+δ1

(2)

φ2=N2+φ2+c(dtr-dts)/λ2+ΔT+ΔI2+δ2

(3)

在式(2)和式(3)中，c是光速，N1和N2分別表示B1和B2的整周模糊度，單位為周。φ表示相位的小數(shù)部分，δ表示其他誤差，包括頻率相關(guān)的接收機和衛(wèi)星端的非校正硬件延遲、相位延遲，觀測噪聲、多路徑效應(yīng)及其他未模型化的誤差。對于同一臺接收機和同一顆衛(wèi)星而言，接收機鐘差dtr、衛(wèi)星鐘差dts，對流層延遲是完全相同的，電離層的誤差由于受到頻率的影響而不同，而且其他誤差δ也具有一定的正相關(guān)性。

本文要組合的波長為λWL=c/(B1-B2)=1.024 7 m，下面對雙頻組合載波相位得到該波長進行推理。設(shè)站星距為ρ，用下標1、2、WL分別表示B1、B2、組合頻率的波長，那么可列出如下關(guān)系

(4)

因為λ1=c/B1，λ2=c/B2，那么式(4)可改寫為

(5)

所以φWL=φ1-φ2，即組合的載波相位可由式(2)減去式(3)得

φWL=(N1+φ1)-(N2+φ2)+ΔI1-ΔI2+

δ1-δ2+c(dtr-dts)(1/λ1-1/λ2)

(6)

設(shè)當(dāng)前歷元為ti，后一歷元為ti+1，若對式(6)中相鄰歷元相減，由于式(6)中的衛(wèi)星鐘差和接收機鐘差基本相同，可以認為基本消去，那么組合的相位觀測量的變化量為

φWL(ti+1)-φWL(ti)=ΔN1-ΔN2+ΔI1-

(7)

把雙頻組合的波長代入式(1)，重構(gòu)出的多普勒為fd。按照最初多普勒積分法探測周跳的思路，結(jié)合式(6)和式(7)，周跳組合可以表示成如下

(8)

式中,ΔN1和ΔN2分別表示B1和B2頻點上發(fā)生的周跳，此時對式(8)右邊得到的結(jié)果進行四舍五入取整，可以得到ΔN1-ΔN2的整數(shù)解。顯然當(dāng)ΔN1=ΔN2時，與無周跳的結(jié)果一樣，因此需要聯(lián)合其他方法把這些周跳盲點也探測出來。

1.3 聯(lián)合探測與修復(fù)周跳

STPIR法是在電離層殘差法上進行時間二次差的方法，當(dāng)電離層活動較大時，能有效減小電離層誤差的影響，探測效果較好[15]。電離層殘差可以表示成

(9)

電離層殘差法是對式(9)進行歷元單差，那么

ΔφPIR(ti)=φPIR(ti+1)-φPIR(ti)

(10)

式中，ti表示當(dāng)前歷元，ti+1表示下一歷元。STPIR法是對式(10)的再次歷元間差分，用ti+2表示當(dāng)前往后第2個歷元，則

φPIR(ti+2)-2φPIR(ti+1)+φPIR(ti)

(11)

在式(10)中，若產(chǎn)生周跳，周跳值應(yīng)為

(12)

因此，如果非連續(xù)歷元，如400歷元產(chǎn)生周跳，那么式(12)的再次差分，即STPIR法就會產(chǎn)生圖2中第400 s、401 s中的大小相等、符號相反的周跳探測量，即分別為ΔN1-f1/f2×ΔN2和f1/f2×ΔN2-ΔN1，此時可取前一個值作為周跳檢測量。但是，若歷元199和歷元200發(fā)生連續(xù)周跳，如圖2中第200 s處(虛線圈出)，就會產(chǎn)生重疊，此時前一個歷元的周跳探測量仍然取前一個值，即2.305周，而后一個歷元的周跳探測量要取后一個值的相反數(shù)，即1.987周作為周跳探測量。

圖2 STPIR法探測周跳Fig.2 STPIR method detect cycle-slip

結(jié)合式(8)和式(12)，可以看出上述兩種方法都存在探測盲點，對式(8)，當(dāng)ΔN1=ΔN2時探測不出；對式(12)，ΔN1≈f1/f2×ΔN2時不敏感，比如(ΔN1,ΔN2)為(9,7)的周跳。那么就需要聯(lián)立方程組解算出周跳。設(shè)重構(gòu)多普勒積分法的周跳探測量為a，STPIR法的周跳探測量為b。由于ΔN1和ΔN2都是整數(shù)，則ΔN1-ΔN2必然是整數(shù)，因此先對a進行四舍五入取整(round)，然后聯(lián)立STPIR法，可得如下方程組

(13)

[ΔN1,ΔN2]T=round(H-1·B)

(14)

周跳的修復(fù)方法：若在周跳探測過程中出現(xiàn)了周跳，把對應(yīng)的相位觀測值φ1和φ2分別減去式(14)得到的跳變值ΔN1和ΔN2即可。

2 重構(gòu)多普勒誤差分析

本節(jié)首先分析用廣播星歷對各類衛(wèi)星重構(gòu)多普勒的影響，然后對比多普勒觀測值和重構(gòu)多普勒的誤差。文獻[16]研究了多普勒觀測值的精度，研究表明多普勒觀測值與衛(wèi)星高度角有關(guān)，高度角越大，誤差越小，可以綜合認為原始多普勒觀測值的精度為2 cm/s。因此，隨著采樣率的減小，即數(shù)據(jù)歷元間隔的增大，多普勒誤差會增大。而分析式(1)可知重構(gòu)多普勒誤差來源于衛(wèi)星的位置、速度，接收機的位置、速度，以及實際的信號波長。

在式(1)中，當(dāng)接收機靜止時v=0，但接收機運動時，需要知道它的運動速度。為了使該方法在動態(tài)中也能進行多普勒估計，推廣其適用性，本文用多普勒觀測值進行接收機速度的估計。在本文所提出方法中加入了用多普勒觀測值計算接收機速度的方法[13,17]，文獻研究表明，大約95%定速結(jié)果可達0.2 m/s的精度。

本文方法已在北斗和GPS衛(wèi)星上進行驗證，GPS只有MEO衛(wèi)星，而北斗導(dǎo)航系統(tǒng)由GEO、MEO、IGSO 3種類型的衛(wèi)星組成[18]，所以此處以北斗為例進行說明。3種類型衛(wèi)星的軌道和速度(相對于地球的速度)都能用廣播星歷的16參數(shù)模型進行計算[19]，北斗2號的MEO和IGSO優(yōu)于2 m[20]，而GEO略差，優(yōu)于5 m[21]，并且北斗3號衛(wèi)星各項性能指標皆優(yōu)于北斗2號[21-22]。按PDOP為5的保守估計，衛(wèi)星速度由16參數(shù)模型計算的誤差應(yīng)小于4 cm/s，而由18參數(shù)模型計算的誤差最多只有10-4m/s量級[23]。目前，偽距單點定位精度可達米級[24]，偽距相對定位精度可達亞米級[25]。

下面用實測靜態(tài)數(shù)據(jù)進行驗證：首先用精密星歷計算得到的衛(wèi)星位置、速度，用精密單點定位確定接收機位置，并進行多普勒重構(gòu)。然后用廣播星歷計算得到衛(wèi)星的位置、速度，用偽距單點定位計算接收機位置，并進行多普勒重構(gòu)。最后以精密星歷重構(gòu)的多普勒為基準，以廣播星歷重構(gòu)的多普勒結(jié)果減去前者。對于廣播星歷得到的衛(wèi)星位置，偽距定位得到的接收機位置，都是米級誤差，這些誤差綜合起來可達十幾米。但是從圖3重構(gòu)多普勒對比的結(jié)果看，GEO衛(wèi)星差別最小，是10-4數(shù)量級，IGSO和MEO是10-3數(shù)量級，這些誤差對重構(gòu)的影響很小。這可以用式(1)中的重構(gòu)原理進行解釋：結(jié)合式(1)和表1，上述米級誤差，經(jīng)過除以107m數(shù)量級的站星距|ρ|，再乘以103m數(shù)量級的衛(wèi)星速度，這些誤差的影響將會降至原來的10-4倍。雖然GEO的軌道精度最差，但是由于它的軌道高度較高，|ρ|較大，且是同步衛(wèi)星，速度相對于其他衛(wèi)星很小，從而使得在重構(gòu)多普勒過程中反而受軌道誤差的影響最小。因此用廣播星歷和偽距單點定位算法得到的衛(wèi)星位置、速度、接收機位置也能較好地重構(gòu)多普勒。而經(jīng)過雙頻組合后的波長變得很長，重構(gòu)其多普勒的誤差會更小。

圖3 各類型衛(wèi)星不同條件下重構(gòu)多普勒的對比Fig.3 Comparison of estimated Doppler for different types of satellites under different conditions

表1 各類衛(wèi)星重構(gòu)多普勒中重要參數(shù)Tab.1 Important parameters in Doppler reconstruction of various satellites

這里再對比北斗B1頻點的多普勒觀測值和重構(gòu)值的周跳探測噪聲：圖4是采用多普勒積分法對一組無周跳的數(shù)據(jù)(30 s歷元間隔)進行周跳檢測的結(jié)果，橫坐標表示觀測歷元，縱坐標表示周跳探測的結(jié)果。理論上無周跳時，探測結(jié)果應(yīng)該小于1，但明顯多普勒觀測值的探測誤差比較大，特別是在圖4(a)中框出的歷元處。對比可知，用重構(gòu)多普勒進行積分探測周跳的效果要優(yōu)于用多普勒觀測值。

圖4 兩種多普勒積分法對比Fig.4 Comparison of two Doppler integral methods

3 試驗與分析

考慮到實際應(yīng)用場景中，動態(tài)采樣率較高，靜態(tài)采樣率較低，本文采用1 Hz采樣率和30 s采樣間隔的實測數(shù)據(jù)分別進行動態(tài)和靜態(tài)模式下的驗證，并且把所提出的重構(gòu)多普勒積分法跟MW法進行對比。為了展示在不同衛(wèi)星類型、不同頻率組合的試驗，其中1 Hz采樣率的數(shù)據(jù)采用GPS的L1和L2頻率組合，30 s采樣間隔的數(shù)據(jù)采用B1和B3頻率組合。

首先用廣播星歷計算出衛(wèi)星位置、衛(wèi)星速度，接著用偽距定位計算出接收機大概位置。選用的數(shù)據(jù)無周跳，在表2的歷元中模擬加入周跳(ΔN1,ΔN2)進行驗證。由于MW的周跳檢測量與重構(gòu)多普勒積分法都是整數(shù)，這里考慮先對其結(jié)果先進行四舍五入取整，再求解，如式(12)、式(13)所示。

本文試驗結(jié)果如圖5—圖10所示，可以看出在無周跳時，重構(gòu)多普勒積分法的探測量比MW法小，即前者探測誤差更小。在動態(tài)模式中，由于引入接收機速度，也引入相關(guān)誤差，會對重構(gòu)多普勒積分法有一定影響，而不會對MW法產(chǎn)生影響。但對比圖5和圖6，顯然還是重構(gòu)多普勒積分法在無周跳時，探測量更接近0，即探測誤差更小。雖然在1 s的高采樣率情況下，MW能準確探測出ΔN1-ΔN2，即兩頻率周跳之差。但在30 s低采樣率時，結(jié)合圖9—圖10以及表2的探測結(jié)果，發(fā)現(xiàn)MW有個別歷元即使無周跳，其探測量也會大于0.5而導(dǎo)致誤判，本文方法在試驗過程中未發(fā)現(xiàn)探測量大于0.5而導(dǎo)致誤判。圖7中框出部分表示相鄰歷元發(fā)生周跳，與圖4中一樣會在中間點產(chǎn)生重疊，采用相同處理方法；試驗結(jié)果表明，可以處理相鄰歷元出現(xiàn)周跳的情況。圖5、圖6中的第450歷元，以及圖9、圖10中的351歷元沒有出現(xiàn)跳變，正說明了相同周跳出現(xiàn)在兩個頻率時，恰好是重構(gòu)多普勒積分法和MW法的共同探測盲點。

表2 兩種聯(lián)合方法的周跳探測結(jié)果

圖5 重構(gòu)多普勒積分法(1 s)Fig.5 Reconstruction Doppler integration method(1 s)

圖6 MW法(1 s)Fig.6 MW method(1 s)

圖7 STPIR法(1 s)Fig.7 STPIR method(1 s)

圖8 STPIR法(30 s)Fig.8 STPIR method(30 s)

圖9 重構(gòu)多普勒積分法(30 s)Fig.9 Reconstruction Doppler integration method(30 s)

圖10 MW法(30 s)Fig.10 MW method(30 s)

文獻[11,26—27]的研究表明，MW法受觀測噪聲(偽距噪聲較大)和多路徑效應(yīng)的影響，導(dǎo)致其探測噪聲較大，都需要對其進行平滑濾波處理才能有好的效果。本文采用的重構(gòu)多普勒的方法誤差主要來源于衛(wèi)星的位置和速度，以及接收機的位置和速度，目前條件下，僅用廣播星歷和偽距定位方法得到這些結(jié)果對該方法影響不大。

圖9中直線外實心點的歷元周跳探測出錯。在歷元302中，MW法的探測結(jié)果是0.707，四舍五入就是1，STPIR法的探測結(jié)果是0.024 9。在歷元375中，MW法的探測結(jié)果是-0.823，四舍五入就是-1，STPIR法的探測結(jié)果是0.005 7。其原因是：B1和B3頻率上的周跳差值都是1，這是由于MW誤判造成的。雖然STPIR法的探測結(jié)果都是正確的，但是根據(jù)式(12)，因為B1/B3=1.230 6，兩頻率上這種比值的周跳是STPIR法的周跳盲點，而5/4=1.25，4/3=1.333，與B1/B3很接近。所以MW的誤判導(dǎo)致用最小二乘法解算得到這些結(jié)果。

4 總 結(jié)

本文提出了一種高精度周跳探測與修復(fù)的方法，與傳統(tǒng)的MW法相比，可更有效地解決周跳探測與修復(fù)的問題。結(jié)果表明，在低采樣率時，重構(gòu)相應(yīng)頻率的多普勒比多普勒觀測值誤差更小。為了更有效地提高重構(gòu)多普勒的周跳檢測精度，進而對雙頻數(shù)據(jù)進行組合，并重構(gòu)其多普勒，進行積分以探測周跳。重構(gòu)多普勒積分法跟MW有共同的探測盲點，而STPIR法也有自身的局限性，用重構(gòu)多普勒積分法聯(lián)合STPIR進行周跳探測。用動態(tài)1 Hz和靜態(tài)1/30 Hz采樣率的實測數(shù)據(jù)，加入模擬周跳進行驗證，并與無平滑濾波的MW法進行對比。結(jié)果表明，重構(gòu)多普勒積分法探測周跳的誤差比MW法小，且在30 s低采樣率時也能準確探測出任意整數(shù)周跳。