劉冰雨, 王中元, 胡 超, 周圣淇

(1.中國礦業(yè)大學(xué) 環(huán)境與測繪學(xué)院,江蘇 徐州 221116; 2.安徽理工大學(xué) 空間信息與測繪工程學(xué)院,安徽 淮南 232001)

2020年6月,最后一顆北斗三號衛(wèi)星成功組網(wǎng),標(biāo)志著北斗衛(wèi)星導(dǎo)航系統(tǒng)(BeiDou Navigation Satellite System,BDS)星座部署全面完成。目前,BDS由8顆地球靜止軌道(Geostationary Earth Orbit,GEO)衛(wèi)星、11顆傾斜地球同步軌道(Inclined Geosynchronous Orbit,IGSO)衛(wèi)星和30顆中地球軌道(Medium Earth Orbit,MEO)衛(wèi)星組成,包含北斗二號全球衛(wèi)星導(dǎo)航系統(tǒng)(BDS-2)和北斗三號全球衛(wèi)星導(dǎo)航系統(tǒng)(BDS-3)2種衛(wèi)星混合星座。相較于BDS-2,BDS-3不僅向下兼容了B1I和B3I信號,還增加了B1C、B2a、B2b和B2(B2a+B2b)信號[1-3]。更多的可用衛(wèi)星和更加豐富的信號資源使得BDS在導(dǎo)航定位中發(fā)揮著更加重要的作用。

差分碼偏差(differential code bias,DCB)是指同一時刻不同頻率或同一頻率上的不同測距信號在發(fā)射鏈路和接收鏈路中所產(chǎn)生時間延遲的差值[4]。DCB在數(shù)量級上可以達(dá)到幾納秒甚至幾十納秒,是影響電離層總電子含量(total electron content,TEC)監(jiān)測和建模的主要誤差源。若忽略DCB的影響,不僅會給TEC的計(jì)算帶來9～30 TECU(TECU為電離層TEC的單位,1 TECU=1×1016個電子/m2)的偏差,也會使得依靠偽距進(jìn)行定位的測量結(jié)果產(chǎn)生數(shù)米偏差[5-6]。常用DCB估計(jì)方法有2種:① 與電離層參數(shù)同步估計(jì)[7-8];② 先用電離層產(chǎn)品直接扣除電離層TEC,再解算DCB[9]。目前,多全球衛(wèi)星導(dǎo)航系統(tǒng)(global navigation satellite system,GNSS)試驗(yàn)跟蹤網(wǎng)(Multi-GNSS Experiment,MGEX)向全球用戶提供2款北斗DCB產(chǎn)品,分別是德國宇航中心(Deutsches Zentrum für Luft-und Raumfahrt,DLR)提供的DLR DCB產(chǎn)品和中國科學(xué)院(Chinese Academy of Sciences,CAS)提供的CAS DCB產(chǎn)品。文獻(xiàn)[10]對MGEX發(fā)布的北斗DCB產(chǎn)品的精度與穩(wěn)定性進(jìn)行質(zhì)量分析,并初步計(jì)算分析了北斗DCB的隨機(jī)誤差特性;文獻(xiàn)[11]分析不同太陽活動水平下BDS衛(wèi)星DCB產(chǎn)品的穩(wěn)定性變化特性,并建模實(shí)現(xiàn)不同太陽活動下BDS衛(wèi)星DCB短期預(yù)報(bào);文獻(xiàn)[12-13]首先明確BDS時間群延遲(timing group delay,TGD)和DCB之間的關(guān)系,然后推導(dǎo)了適用多種場合的北斗DCB改正模型,并論證在高精度定位中采用合適DCB改正模型的重要性;文獻(xiàn)[14]在文獻(xiàn)[12]的基礎(chǔ)上推導(dǎo)三頻無電離層組合DCB改正模型,并用實(shí)驗(yàn)分析驗(yàn)證DCB改正對標(biāo)準(zhǔn)單點(diǎn)定位(standard point positioning,SPP)和精密單點(diǎn)定位2種定位模式的影響。

以上關(guān)于北斗DCB穩(wěn)定性變化特性和DCB改正定位模型的分析,都基于DCB為1 d中的常量參數(shù),關(guān)于BDS衛(wèi)星DCB短時變化特性的研究很少。本文分析BDS衛(wèi)星DCB在1 d中的短時變化,針對單歷元數(shù)據(jù)進(jìn)行短時DCB估計(jì),即在每個觀測歷元估計(jì)出該時刻對應(yīng)的DCB估值。首先采用最小二乘和Tikhonov正則化方法,解算BDS衛(wèi)星在各個歷元時刻下的DCB估值,然后對得到的各衛(wèi)星DCB序列的精度、穩(wěn)定性進(jìn)行分析,并采用譜分析方法確定短時DCB周期性,構(gòu)建模型擬合短時DCB變化,最后通過實(shí)驗(yàn)來探討短時DCB序列改正對SPP的影響。

1 短時DCB估計(jì)方法

通過雙頻無幾何組合,可得到偽距和載波相位無幾何觀測量,表達(dá)式為:

(1)

(2)

由于偽距無幾何觀測量觀測噪聲較大,常用載波相位對其進(jìn)行平滑處理,平滑后的偽距觀測值可表示為:

(3)

其中:P4,sm為平滑后的偽距無幾何觀測量;P4,prd(t)為偽距和載波的組合觀測量;ωt為與第t個歷元相關(guān)的比例系數(shù),ωt=1/t。當(dāng)t=1時,P4,sm=P4。

電離層延遲通常忽略高階項(xiàng)影響,采用一階表達(dá)式表示,即

(4)

其中:fj為信號頻率,j=1,2;Stec為傾斜TEC。將(4)式代入(3)式,可得:

(5)

由(5)式可得Stec的表達(dá)式為:

(6)

為方便電離層建模,常借助一個映射函數(shù)FM將傾斜方向Stec轉(zhuǎn)到天頂方向Vtec,表達(dá)式為:

(7)

其中:Vtec為天頂方向TEC;z為衛(wèi)星高度角的余角;R為地球平均半徑;H為電離層薄層高度,本文取值為506.7 km;α為比例因子,本文取值為0.978 2。將(7)式代入(6)式,可得:

(P4,sm-cDs-cDr)

(8)

本文采用球諧函數(shù)模型對天頂方向Vtec進(jìn)行建模,表達(dá)式為:

[anmcos(ms)+bnmsin(ms)]

(9)

將(8)式、(9)式聯(lián)立可得:

(10)

(10)式改寫成矩陣形式為:

(11)

在解算短時DCB過程中,鄰近測站對應(yīng)穿刺點(diǎn)緯度和日固經(jīng)度變化較小,系數(shù)矩陣中數(shù)據(jù)存在一定的相關(guān)性,在某些歷元中直接采用最小二乘法估算會出現(xiàn)秩虧現(xiàn)象。為了減弱或消除病態(tài)問題對秩虧歷元的影響,本文采用Tikhonov正則化[15]對秩虧歷元的參數(shù)進(jìn)行解算。Tikhonov正則化方法是對最小二乘法的改進(jìn),正則化最小二乘代價(jià)函數(shù)為:

(12)

(12)式的解為:

(13)

其中,λ為正則化參數(shù),λ≥0,當(dāng)λ=0時該解為最小二乘解。為了選取合適的λ,本文采用嶺跡法確定秩虧歷元中λ參數(shù),對隨機(jī)選取的50個秩虧歷元繪制嶺跡圖,當(dāng)λ處于0.001附近時,嶺跡趨于平穩(wěn),故本文解算過程中λ取0.001。

2 數(shù)據(jù)處理與實(shí)驗(yàn)分析

本文選取2021年年積日第152天至第212天國際GNSS服務(wù)(International GNSS Service,IGS)和MGEX測站網(wǎng)BDS的C2I、C6I觀測數(shù)據(jù),測站分布位置如圖1所示。采用4階球諧函數(shù)模型對TEC進(jìn)行建模,衛(wèi)星截止高度角設(shè)置為10°。為了分析DCB短時變化,同步估計(jì)每個觀測歷元時刻下的DCB參數(shù)和球諧函數(shù)參數(shù)。針對可能出現(xiàn)的極少數(shù)歷元缺失情況,采用前后2個歷元平均值進(jìn)行替代。

圖1 測站位置分布圖

以年積日第200天為例,按照本文方法解算出的部分衛(wèi)星DCB估值序列如圖2所示。從圖2可以看出:BDS-2中C02、C04、C06、C08衛(wèi)星在各個歷元的估值比較穩(wěn)定,僅有輕微的波動,1 d中上下波動在0.5 ns左右;BDS-3中C26、C32、C34、C40衛(wèi)星在某些時段會出現(xiàn)較大的起伏,穩(wěn)定性與BDS-2中的衛(wèi)星相比較差,1 d中上下波動在1.2 ns左右。

圖2 部分衛(wèi)星短時DCB序列

2.1 BDS衛(wèi)星短時DCB精度與穩(wěn)定性分析

為了驗(yàn)證BDS衛(wèi)星短時DCB估計(jì)精度,將CAS發(fā)布的DCB產(chǎn)品作為參考值,以本文方法得到的1 d中各衛(wèi)星短時DCB序列分別與相應(yīng)的參考值作差,可獲得各衛(wèi)星DCB估值與參考值的差值序列,對各衛(wèi)星的差值序列求取平均值,可得到各衛(wèi)星短時DCB序列在1 d中的精度狀況?？紤]到僅對比1 d數(shù)據(jù)具有偶然性,對2021年年積日第152天至第212天共60 d天數(shù)據(jù)進(jìn)行絕對值平均。目前,BDS-2衛(wèi)星有15顆,偽隨機(jī)噪聲(pseudo-random noise,PRN)碼編號為C01～C16(C15編號未用);BDS-3衛(wèi)星PRN碼編號為C19～C46,其中C31衛(wèi)星處于試驗(yàn)階段,未參與解算。各衛(wèi)星在60 d中的平均偏差如圖3所示。

圖3 BDS衛(wèi)星短時DCB(C2I、C6I)平均偏差

從圖3可以看出,BDS-2衛(wèi)星DCB短時估計(jì)的精度小于BDS-3衛(wèi)星。BDS-2衛(wèi)星短時DCB估值與參考值在60 d中的平均偏差為0.34 ns,且各衛(wèi)星平均偏差之間相差較大,C04衛(wèi)星與參考值的平均偏差(0.79 ns)最大,C16衛(wèi)星與參考值的平均偏差(0.11 ns)最小。BDS-3衛(wèi)星短時DCB在60 d 中的平均偏差為0.14 ns,各衛(wèi)星短時DCB平均偏差相差較小,C38衛(wèi)星平均偏差(0.35 ns)最大, C34衛(wèi)星平均偏差(0.08 ns)最小,其他衛(wèi)星平均偏差大多在0.10 ns左右。

為了分析衛(wèi)星短時DCB在日內(nèi)變化的穩(wěn)定性,本文統(tǒng)計(jì)各觀測日1 d內(nèi)衛(wèi)星短時DCB的標(biāo)準(zhǔn)差(standard deviation,STD)σ,并將其作為反映衛(wèi)星短時DCB穩(wěn)定性的指標(biāo)[9]。σ表達(dá)式為:

(14)

圖4 BDS衛(wèi)星平均日內(nèi)變化穩(wěn)定性

BDS-2中GEO衛(wèi)星包括C01～C05衛(wèi)星,BDS-2中非GEO衛(wèi)星包括C06～C16衛(wèi)星(C15編號未用),BDS-3中非GEO衛(wèi)星為C19～C46衛(wèi)星。從圖4可以看出:BDS-2和BDS-3中的非GEO衛(wèi)星(MEO衛(wèi)星和IGSO衛(wèi)星)的平均日內(nèi)穩(wěn)定性非常接近,其STD均值分別為0.42、0.43 ns；BDS-2中的GEO衛(wèi)星日內(nèi)變化穩(wěn)定性更高,其STD均值為0.24 ns。

2.2 短時DCB周期性分析與擬合建模

通過對年積日第152天至第212天各衛(wèi)星DCB短時序列進(jìn)行譜分析,可發(fā)現(xiàn)DCB短時估值的周期性變化。本文使用快速傅里葉變換對DCB序列進(jìn)行分析,得到部分衛(wèi)星在各觀測日中的頻譜圖,如圖5所示。

從圖5a可以看出,BDS的GEO衛(wèi)星DCB估值在頻率為2.035×10-6、1.170×10-5、2.289×10-5Hz會達(dá)到一個極大值,此時三者對應(yīng)的周期分別為5.6 d、23.7 h、12.1 h,即GEO衛(wèi)星會表現(xiàn)出以5.6 d、23.7 h、12.1 h為周期的周期性變化。

從圖5b可以看出:IGSO衛(wèi)星在頻率為1.157×10-5、2.314×10-5、3.471×10-5Hz或4.641×10-5Hz會達(dá)到一個極大值；在頻率為2.314×10-5Hz時振幅都會達(dá)到最大,在1.157×10-5Hz時振幅次之,上述2個頻率對應(yīng)的周期分別為12.0、24.0 h;個別衛(wèi)星在頻率為3.471×10-5Hz 或4.641×10-5Hz時也會達(dá)到一個相對較小的峰值,上述2個頻率對應(yīng)的周期分別為8.0、6.0 h。這說明IGSO衛(wèi)星除了會有以12.0、24.0 h為周期的強(qiáng)周期性變化,也會在某衛(wèi)星上表現(xiàn)出以6.0、8.0 h為周期的弱周期性變化。

從圖5c可以看出:MEO衛(wèi)星在頻率為9.918×10-6、2.149×10-5、3.319×10-5Hz或4.311×10-5Hz會達(dá)到一個極大值；在頻率為9.918×10-6、4.311×10-5Hz時振幅較為明顯,上述2個頻率對應(yīng)的周期分別為28.0、6.4 h;個別衛(wèi)星在頻率為2.149×10-5Hz或3.319×10-5Hz時也會達(dá)到一個相對較小的峰值,此時對應(yīng)的周期分別為12.9、8.4 h。因此,MEO衛(wèi)星有以6.4、28.0 h為周期的強(qiáng)周期性變化,也會出現(xiàn)以12.9、8.4 h為周期的弱周期性變化。

圖5 BDS部分衛(wèi)星短時DCB序列頻譜

由上述周期性分析可知,衛(wèi)星短時DCB在1 d內(nèi)存在周期性變化,考慮到此特性,可將短時DCB的擬合模型表示為一個二次函數(shù)與幾個周期函數(shù)的總和。

表達(dá)式為:

y=A0T2+B0T+C0+

(15)

其中:y為DCB值;T為相對起始?xì)v元的時間間隔;A0、B0、C0分別為二次函數(shù)的二次項(xiàng)系數(shù)、一次項(xiàng)系數(shù)和常數(shù)項(xiàng);Dk、Ek為周期函數(shù)在第k個周期的系數(shù);Qk為周期函數(shù)的周期。為了得到每日短時DCB的擬合函數(shù),首先對每日衛(wèi)星DCB短時估值進(jìn)行譜分析,確定周期函數(shù)的周期后,采用最小二乘法確定周期函數(shù)和二次函數(shù)的系數(shù)。

以年積日第200天為例,部分衛(wèi)星DCB擬合函數(shù)中系數(shù)與周期函數(shù)周期分別見表1、表2所列。

表1 擬合模型中二次函數(shù)對應(yīng)系數(shù)

表2 擬合模型中周期函數(shù)對應(yīng)系數(shù)與周期

為了評估模型的擬合精度,本文將擬合模型得到的擬合值與解算出的DCB估值進(jìn)行比較。部分衛(wèi)星短時DCB采用(15)式的擬合值與解算值對比如圖6所示。擬合值與解算值差值的相關(guān)統(tǒng)計(jì)特性見表3所列。表3中,RMS表示差值的均方根(root mean square)值。

從圖6可以看出,各衛(wèi)星DCB擬合值與DCB解算值具有很好的一致性。

從表3可以看出,各衛(wèi)星DCB擬合值與解算值的平均差值幾乎為0 ns,所選6顆衛(wèi)星的RMS值分別為0.006、0.017、0.053、0.006、0.018、0.047 ns,表明本文提出的擬合模型具有較高的精度,可以很好地?cái)M合衛(wèi)星DCB在日內(nèi)的短時變化,也從側(cè)面反映了采用譜分析方法分析周期性的可行性。

圖6 短時衛(wèi)星DCB擬合結(jié)果

表3 解算值與擬合值之間差值統(tǒng)計(jì)特性單位:ns

3 DCB改正對SPP的影響

為了分析短時DCB估計(jì)對SPP的影響,選用IGS測站網(wǎng)下ABPO、ULAB、LEIJ和WUH2測站在2021年年積日第200天的數(shù)據(jù)進(jìn)行SPP實(shí)驗(yàn)。對BDS的B1頻點(diǎn)觀測量分別采用3種不同的處理策略,即不改正DCB、CAS產(chǎn)品改正、本文方法改正,分析單頻SPP結(jié)果。DCB改正的單頻SPP模型見文獻(xiàn)[12]。4個測站在3種策略下單頻SPP三維(3D)殘差結(jié)果如圖7所示。

從圖7可以看出:

(1) 經(jīng)過DCB改正后,單頻SPP的定位精度提升明顯,基本在5 m以內(nèi)。

(2) 使用CAS產(chǎn)品與本文方法對單頻SPP改正,3D殘差的量級和歷元間變化趨勢大致相當(dāng)。在不同測站上,兩者定位精度互有優(yōu)劣,互差在cm級,改正效果大致相當(dāng)。

4個測站單頻SPP在不同方向上的定位結(jié)果見表4所列。

由表4可知:

(1) CAS產(chǎn)品改正和本文方法改正對SPP的影響為m級,CAS產(chǎn)品改正的平面方向、U方向和3D方向平均RMS分別為1.118、3.162、3.407 m,本文方法改正的平面方向、U方向和3D方向平均RMS分別為1.112、3.143、3.387 m?？梢钥闯?兩者改正后SPP定位精度互差在cm級,在平面方向、U方向和3D方向,本文方法改正與CAS產(chǎn)品改正相比,略微提升0.6、1.9、2.0 cm,其原因可能是短時DCB參數(shù)更新周期短,能更好反映DCB在各個歷元的真實(shí)狀況。

(2) CAS產(chǎn)品改正的SPP在E、N、U和3D方向定位精度平均提升55.5%、58.6%、42.5%、45.0%；本文方法改正的SPP在E、N、U和3D方向定位精度平均提升56.9%、58.2%、43.0%、45.3%。

圖7 4個測站DCB改正前、后SPP定位3D殘差對比

表4 4個測站3種DCB處理策略下SPP定位平均RMS值

4 結(jié) 論

通過在單歷元間估算DCB,得到衛(wèi)星DCB短時序列,可以對DCB短時序列變化特性進(jìn)行分析,也能充分考慮DCB短時變化改正對導(dǎo)航定位的影響。

本文在各個歷元間采用最小二乘和Tikhonov正則化方法同步估算電離層參數(shù)和BDS的DCB短時序列,分析衛(wèi)星DCB短時精度、穩(wěn)定性和周期性,最后通過SPP實(shí)驗(yàn)分析DCB短時序列對定位的影響。結(jié)果表明：施加DCB短時改正和施加CAS DCB產(chǎn)品都能顯著提高SPP定位精度,定位精度提升都在40%以上;施加短時DCB改正和CAS產(chǎn)品改正,在不同測站上SPP定位精度各有優(yōu)劣,兩者差異在cm級,改正效果大致相當(dāng)。