多普勒積分重構(gòu)與STPIR聯(lián)合周跳探測(cè)與修復(fù)

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

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

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

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

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

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

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

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

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

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

(1)

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

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

分析式(1)可知重構(gòu)多普勒誤差主要來(lái)源于中分子計(jì)算的部分，若組合出來(lái)的波長(zhǎng)越長(zhǎng)，即式(1)中分母變大，重構(gòu)誤差就會(huì)減小。下面以北斗為例，北斗的頻率B1=1 561.098 MHz和B2=1 207.14 MHz，對(duì)應(yīng)的相位觀測(cè)值φ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)的接收機(jī)和衛(wèi)星端的非校正硬件延遲、相位延遲，觀測(cè)噪聲、多路徑效應(yīng)及其他未模型化的誤差。對(duì)于同一臺(tái)接收機(jī)和同一顆衛(wèi)星而言，接收機(jī)鐘差dtr、衛(wèi)星鐘差dts，對(duì)流層延遲是完全相同的，電離層的誤差由于受到頻率的影響而不同，而且其他誤差δ也具有一定的正相關(guān)性。

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

(4)

因?yàn)棣?=c/B1，λ2=c/B2，那么式(4)可改寫(xiě)為

(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，若對(duì)式(6)中相鄰歷元相減，由于式(6)中的衛(wèi)星鐘差和接收機(jī)鐘差基本相同，可以認(rèn)為基本消去，那么組合的相位觀測(cè)量的變化量為

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

(7)

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

(8)

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

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

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

(9)

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

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

(10)

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

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

(11)

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

(12)

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

圖2 STPIR法探測(cè)周跳Fig.2 STPIR method detect cycle-slip

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

(13)

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

(14)

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

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

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

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

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

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

圖3 各類型衛(wèi)星不同條件下重構(gòu)多普勒的對(duì)比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

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

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

3 試驗(yàn)與分析

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

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

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

表2 兩種聯(lián)合方法的周跳探測(cè)結(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)

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

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

4 總 結(jié)

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