于 燁,黃 默,段 濤,王長(zhǎng)元,胡 銳
(1.中國(guó)科學(xué)院 微電子研究所,北京 100029; 2.中國(guó)科學(xué)院大學(xué),北京 100049;3.中國(guó)科學(xué)院大學(xué) 微電子學(xué)院,北京 100049)
衛(wèi)星鐘差預(yù)報(bào)是導(dǎo)航定位中一項(xiàng)至關(guān)重要的工作之一[1-2].星載原子鐘因其自身的性能以及所處空間復(fù)雜環(huán)境的影響,鐘源偏差需要與地面站的原子鐘進(jìn)行比對(duì)并校準(zhǔn)以保持比較精確的時(shí)間信息,但是出于某些原因星載原子鐘無(wú)法與地面站保持比對(duì)的時(shí)候就需要采用已有的鐘差信息進(jìn)行預(yù)報(bào).在這段時(shí)間內(nèi),衛(wèi)星以預(yù)報(bào)的星歷信息進(jìn)行自主完成軌道確定和廣播星歷的播發(fā),這對(duì)于衛(wèi)星的自主生存能力具有重要影響[3-4],所以如何提高衛(wèi)星鐘差預(yù)報(bào)的精度和穩(wěn)定度一直以來(lái)都是國(guó)內(nèi)外研究的熱點(diǎn)問(wèn)題之一.
長(zhǎng)久以來(lái),在衛(wèi)星鐘差預(yù)報(bào)方面,國(guó)內(nèi)外許多學(xué)者進(jìn)行了多角度、多方位和深層次的研究,并取得了一系列重要研究成果[5-7].其中,灰色系統(tǒng)預(yù)報(bào)模型只需要少量的鐘差數(shù)據(jù)建模就可以預(yù)報(bào),并且需要數(shù)據(jù)平滑且呈類(lèi)指數(shù)變化,對(duì)于數(shù)據(jù)呈非指數(shù)變化的鐘差預(yù)報(bào)效果較差等特點(diǎn).文獻(xiàn)[8]提出采用最小一乘法準(zhǔn)則去估計(jì)灰色模型的參數(shù),以克服在鐘差波動(dòng)較大的情況下最小二乘法準(zhǔn)則估計(jì)參數(shù)預(yù)報(bào)精度的不足,并且通過(guò)預(yù)報(bào)試驗(yàn)證明了算法的有效性.文獻(xiàn)[9]提出了一種先用對(duì)數(shù)平滑處理的方法對(duì)原始鐘差數(shù)據(jù)序列進(jìn)行平滑處理,再使用自適應(yīng)雙子群改進(jìn)粒子群算法來(lái)優(yōu)化灰色模型的發(fā)展系數(shù)和灰作用量的預(yù)報(bào)策略,從而增強(qiáng)了灰色模型的泛化能力.文獻(xiàn)[10]提出先對(duì)衛(wèi)星鐘差第一天偶數(shù)整時(shí)的值采用灰色模型進(jìn)行建模,然后利用已建立的模型預(yù)報(bào)第二天相應(yīng)時(shí)間的鐘差,對(duì)預(yù)報(bào)得到的鐘差進(jìn)行拉格朗日插值建模,通過(guò)內(nèi)插得到其他時(shí)刻的預(yù)報(bào)值,以此來(lái)提高模型的預(yù)報(bào)能力.文獻(xiàn)[11]分析了灰色模型發(fā)展系數(shù)和灰作用量對(duì)預(yù)報(bào)精度的影響,之后通過(guò)微調(diào)這兩個(gè)參數(shù)來(lái)提高灰色模型的預(yù)報(bào)性能,但是針對(duì)鐘差數(shù)據(jù)序列呈現(xiàn)出不同變化趨勢(shì)時(shí),沒(méi)有給出有效的參數(shù)優(yōu)化方法.
從以上的分析可知,研究者們主要是在優(yōu)化灰色模型的兩個(gè)參數(shù)方面對(duì)模型進(jìn)行改進(jìn),沒(méi)有從灰色模型的一次累加生成序列的預(yù)報(bào)模型方面進(jìn)行其他處理.基于此,本文提出了一種基于粒子群算法優(yōu)化指數(shù)函數(shù)和線性函數(shù)逼近加權(quán)灰色回歸組合的自適應(yīng)衛(wèi)星鐘差預(yù)報(bào)方法.在建模之前,考慮到衛(wèi)星鐘差鐘跳頻繁的現(xiàn)象,首先采用中位數(shù)法對(duì)衛(wèi)星鐘差進(jìn)行異常鐘差探測(cè),將異常鐘差剔除后采用分段線性插值法將缺失的鐘差補(bǔ)齊.然后考慮到衛(wèi)星鐘差存在系統(tǒng)噪聲,采用三點(diǎn)平滑法對(duì)鐘差數(shù)據(jù)進(jìn)行平滑處理后建立了粒子群算法優(yōu)化指數(shù)函數(shù)和線性函數(shù)逼近加權(quán)灰色回歸組合的自適應(yīng)衛(wèi)星鐘差預(yù)報(bào)模型,并結(jié)合4種典型變化趨勢(shì)的衛(wèi)星鐘差序列,采用IGS服務(wù)器上發(fā)布的事后精密衛(wèi)星鐘差產(chǎn)品作為試驗(yàn)數(shù)據(jù),進(jìn)行了預(yù)報(bào)試驗(yàn).
由于衛(wèi)星鐘差數(shù)據(jù)跳變頻繁,跳變的鐘差數(shù)據(jù)很不利于建模,所以在建模之前對(duì)衛(wèi)星鐘差數(shù)據(jù)的質(zhì)量進(jìn)行檢測(cè)就顯得很重要.本文采用時(shí)效性和抗差性較好的中位數(shù)法(Median Absolute Deviation, MAD)進(jìn)行異常鐘差探測(cè)[12],其表達(dá)式為
(1)
式中Median(yi)為衛(wèi)星鐘差頻率數(shù)據(jù)yi的中間數(shù).
當(dāng)衛(wèi)星鐘差的頻率yi>(m+n·MAD)或yi<(m+n·MAD)時(shí),可以判斷該點(diǎn)為異常鐘差.將探測(cè)到的異常鐘差剔除后,采用分段線性插值法把缺失的鐘差數(shù)據(jù)補(bǔ)齊.
由于星載原子鐘受自身復(fù)雜的時(shí)頻特性和外界復(fù)雜環(huán)境的影響,穩(wěn)定性存在較大差異的特點(diǎn),而衛(wèi)星鐘差是否平滑對(duì)于模型的預(yù)報(bào)性能有很大影響,所以對(duì)觀測(cè)的鐘差進(jìn)行平滑性檢測(cè)就顯得很有必要.
若ρ(i)滿(mǎn)足:
(2)
則稱(chēng)鐘差數(shù)據(jù)序列X(1)為平滑序列;反之,稱(chēng)鐘差數(shù)據(jù)序列X(1)為非平滑序列.
由于衛(wèi)星鐘差數(shù)據(jù)存在系統(tǒng)噪聲,采用三點(diǎn)平滑法處理,可以提高衛(wèi)星鐘差數(shù)據(jù)序列的平滑度.三點(diǎn)平滑法是一種取算術(shù)平均的方法,它通過(guò)重新分配待處理的衛(wèi)星鐘差與前后鐘差數(shù)據(jù)的權(quán)值,以加強(qiáng)待處理鐘差數(shù)據(jù)的權(quán)重,從而減小衛(wèi)星鐘差的波動(dòng)性,增強(qiáng)待處理的衛(wèi)星鐘差與前后鐘差數(shù)據(jù)間的聯(lián)系.其過(guò)程如下:
設(shè)有一組衛(wèi)星鐘差數(shù)據(jù)序列為
X(1)={x(1)(1),x(1)(2),…,x(1)(n)},
(3)
平滑處理位于衛(wèi)星鐘差中間的數(shù)據(jù)為
(4)
而對(duì)于鐘差數(shù)據(jù)序列兩端的數(shù)據(jù)平滑處理的計(jì)算為
(5)
經(jīng)式(4)和(5)處理后,得到一組噪聲減小的鐘差數(shù)據(jù)序列,然后建模預(yù)報(bào).
灰色模型是指信息不完全知道的系統(tǒng),具體是指它的部分信息是已知的、其余的信息為未知的一種系統(tǒng).其最常見(jiàn)的灰色系統(tǒng)是GM(1,1)模型,它由一個(gè)包含單變量的一階微分方程所構(gòu)成,可以實(shí)現(xiàn)對(duì)自身數(shù)據(jù)的預(yù)測(cè)[13].
設(shè)X(0)={x(0)(1),x(0)(2),…,x(0)(n)}為不同時(shí)刻的衛(wèi)星鐘差,經(jīng)過(guò)中位數(shù)法探測(cè)剔除異常鐘差數(shù)據(jù),然后采用分段線性插值法把缺失的鐘差數(shù)據(jù)補(bǔ)齊,得到一組新的鐘差數(shù)據(jù)序列為X(1)={x(1)(1),x(1)(2),…,x(1)(n)}.再使用三點(diǎn)平滑法處理,以得到一組噪聲減小的鐘差數(shù)據(jù)序列為X(2)={x(2)(1),x(2)(2),…,x(2)(n)},然后對(duì)X(2)這組鐘差數(shù)據(jù)序列建立灰色模型,則這組數(shù)據(jù)的一次累加鐘差數(shù)據(jù)序列為X(3)={x(3)(1),x(3)(2),…,x(3)(n)}的灰色預(yù)報(bào)模型可表示為
(6)
記參數(shù)為A=[au]T,根據(jù)最小二乘法可得
(7)
式中
將式(7)代入式(6)即可預(yù)報(bào)出鐘差的一次累加序列,然后還原即可得到原始鐘差的預(yù)報(bào)值.
分析可以看出式(6)為以下形式
(8)
由于上式為指數(shù)函數(shù)和線性函數(shù)的疊加形式,下面采用指數(shù)曲線y=aex與線性曲線y=αx+β這兩種模型的和來(lái)逼近鐘差序列X(2)的一次累加生成的序列X(3)的預(yù)報(bào)值,如下
(9)
式中,參數(shù)v和參數(shù)C1、C2、C3為待求的參數(shù).
為確定式(9)中的4個(gè)未知參數(shù),構(gòu)造如下序列:
(10)
設(shè):
(11)
則:
(12)
用式(12)除以式(11)得:
(13)
將式(13)兩邊取對(duì)數(shù)得:
(14)
當(dāng)m=1時(shí):
(15)
當(dāng)m=2時(shí):
(16)
直到當(dāng)m=n-3時(shí):
(17)
(18)
(19)
由于離預(yù)報(bào)值近一點(diǎn)的鐘差對(duì)預(yù)報(bào)的鐘差影響大一點(diǎn),離預(yù)報(bào)值遠(yuǎn)一點(diǎn)的鐘差對(duì)預(yù)報(bào)的鐘差影響小一點(diǎn),下面采用變權(quán)加權(quán)法來(lái)解式(19).變權(quán)加權(quán)法相對(duì)于普通加權(quán)法最大的不同是,變權(quán)加權(quán)法對(duì)鐘差時(shí)間序列的可靠性成比例的分別賦予不同的權(quán)值[14],具體加權(quán)方法如下:
Pk=Ri-1,i=1,2,…,n.
(20)
式中:R稱(chēng)為精度遞增因子且R∈[1,2].
綜合式(19)和式(20),利用最小二乘法可得
(21)
式中
精度遞增因子R的選取是該算法預(yù)報(bào)性能好壞的關(guān)鍵,本文采用粒子群優(yōu)化算法(Particle Swarm Optimization,PSO)對(duì)精度遞增因子R進(jìn)行自適應(yīng)尋優(yōu).
PSO算法是由美國(guó)兩位科學(xué)家受鳥(niǎo)類(lèi)群體尋找食物的行為得到啟發(fā)而提出的一種智能優(yōu)化算法.因?yàn)镻SO算法收斂速度比較快,調(diào)整參數(shù)簡(jiǎn)單,所以得到了廣泛的應(yīng)用.因此,本文選用PSO算法對(duì)上述精度遞增因子R進(jìn)行自適應(yīng)尋優(yōu),PSO算法原理[15]如下:
假設(shè)在一個(gè)目標(biāo)搜索空間,其維數(shù)為D維且有N個(gè)粒子組成一個(gè)群落,其中第i個(gè)粒子可表示為一個(gè)D維向量
Xi=(xi1,xi2,…,xiD),i=1,2,…,N.
(22)
第i個(gè)粒子的“飛行”速度也是一個(gè)D維向量,記為
Vi=(vi1,vi2,…,viD),i=1,2,…,N.
(23)
第i個(gè)粒子迄今為止搜索到的最優(yōu)位置稱(chēng)為個(gè)體極值,記為
pbest=(pi1,pi2,…,piD),i=1,2,…,N.
(24)
整個(gè)粒子群迄今為止搜索到的最優(yōu)位置為全局最優(yōu)值,記為
gbest=(pg1,pg2,…,pgD).
(25)
在未找到這兩個(gè)最優(yōu)值時(shí),粒子根據(jù)如下的公式來(lái)更新自己的速度和位置:
(26)
式中:c1和c2為加速常數(shù),一般取c1=c2=2;r1和r2為[0,1]范圍內(nèi)的均勻隨機(jī)數(shù);wi為慣性權(quán)重,wmax和wmin為預(yù)先確定的最大和最小慣性權(quán)重,一般wmax取0.9,wmin取0.1;Gmax為最大迭代次數(shù);Gi為當(dāng)前迭代次數(shù).
采用PSO算法對(duì)精度遞增因子R進(jìn)行自適應(yīng)尋優(yōu)的具體步驟如下:
1)初始化粒子群的速度、位置、種群規(guī)模、最大迭代次數(shù)、最大最小慣性權(quán)重、訓(xùn)練樣本數(shù)m和加速常數(shù),同時(shí)定義粒子群的適應(yīng)度函數(shù)為
(27)
式中:yi為預(yù)報(bào)的鐘差,yip為實(shí)際的觀測(cè)鐘差.
2)計(jì)算各粒子的適應(yīng)度函數(shù)值fi,然后再比較它們的適應(yīng)度函數(shù)值fi的大小.如果粒子的當(dāng)前適應(yīng)度值fi小于該粒子的歷史最優(yōu)適應(yīng)度值fibest,則把fibest賦給fi;同理若fibest小于當(dāng)前粒子種群的最優(yōu)適應(yīng)度值fgbest,則把fgbest賦給fibest.
3)根據(jù)式(26)更新粒子的速度、位置和慣性權(quán)重,然后判斷是否達(dá)到了最大迭代次數(shù)Gmax.若達(dá)到,則迭代結(jié)束,所得到的最優(yōu)位置即為精度遞增因子R的值;若未達(dá)到,則繼續(xù)迭代直到迭代結(jié)束.此算法的具體流程見(jiàn)圖1.
圖1 粒子群優(yōu)化加權(quán)灰色回歸組合自適應(yīng)預(yù)報(bào)算法流程
為了驗(yàn)證本文預(yù)報(bào)方法的有效性和可行性,采用IGS服務(wù)器(ftp://igs.ign.fr/pub/igs/)上發(fā)布的GPS 2 000周,間隔時(shí)長(zhǎng)為15 min的SP3格式的事后精密衛(wèi)星鐘差數(shù)據(jù)進(jìn)行預(yù)報(bào)試驗(yàn).在該時(shí)間段內(nèi)在軌的GPS衛(wèi)星有31顆,其星載原子鐘有BLOCK IIR-Rb鐘、BLOCK IIR-M-Rb鐘、BLOCK IIF-Rb鐘和BLOCK IIF-Cs四種類(lèi)型.由于北斗衛(wèi)星導(dǎo)航系統(tǒng)的星載鐘與GPS基本一致,且北斗二代衛(wèi)星導(dǎo)航系統(tǒng)均搭載的為銣原子鐘,為使研究結(jié)果能為我國(guó)北斗衛(wèi)星導(dǎo)航系統(tǒng)在鐘差預(yù)報(bào)方面提供一些參考,所以選取了GPS IIF-Rb PRN03、GPS IIR-M-Rb PRN12、GPS IIR-Rb PRN14和GPS IIR-M-Rb PRN17四顆搭載銣原子鐘的衛(wèi)星的鐘差數(shù)據(jù)進(jìn)行預(yù)報(bào)試驗(yàn).截止2019年10月16日,它們的相關(guān)信息見(jiàn)表1.
表1 選擇的衛(wèi)星相關(guān)信息
這四顆衛(wèi)星第2 000周第0 d的事后精密鐘差時(shí)間序列的變化情況如圖2所示,其中PRN03號(hào)衛(wèi)星的鐘差時(shí)間序列呈正值單調(diào)遞增趨勢(shì),PRN12號(hào)衛(wèi)星的鐘差時(shí)間序列呈正值單調(diào)遞減趨勢(shì),PRN14號(hào)衛(wèi)星的鐘差時(shí)間序列呈負(fù)值遞減趨勢(shì),PRN17號(hào)衛(wèi)星的鐘差時(shí)間序列呈負(fù)值單調(diào)遞增趨勢(shì),具有充分的代表性.
圖2 原始觀測(cè)衛(wèi)星鐘差變化趨勢(shì)
由于衛(wèi)星鐘差的有效位數(shù)比較多,使得原始觀測(cè)鐘差數(shù)據(jù)異常點(diǎn)容易被掩蓋,而異常鐘差在其對(duì)應(yīng)的頻率數(shù)據(jù)中表現(xiàn)為顯著的峰值點(diǎn),所以先將衛(wèi)星鐘差轉(zhuǎn)換為相應(yīng)的頻率數(shù)據(jù)后更容易對(duì)異常鐘差進(jìn)行探測(cè),具體轉(zhuǎn)換公式如下
(28)
式中:yi表示衛(wèi)星鐘差的頻率數(shù)據(jù),xi+1和xi分別表 示第i+1和i歷元的衛(wèi)星鐘差值(相位數(shù)據(jù)),τ0表示采樣間隔.
PRN03、PRN12、PRN14和PRN17號(hào)衛(wèi)星相應(yīng)的頻率數(shù)據(jù)變化見(jiàn)圖3,經(jīng)中位數(shù)法探測(cè)發(fā)現(xiàn),在此時(shí)間段內(nèi),PRN03、PRN14和PRN17三顆衛(wèi)星存在異常的鐘差數(shù)據(jù),將探測(cè)的異常鐘差數(shù)據(jù)剔除后采用分段線性插值法將缺失的數(shù)據(jù)補(bǔ)齊,補(bǔ)齊后的變化見(jiàn)圖4,然后還原為鐘差的相位數(shù)據(jù)進(jìn)行建模預(yù)報(bào).
IGS發(fā)布的超快速衛(wèi)星星歷產(chǎn)品的更新時(shí)間約為6 h,鑒于此,為了驗(yàn)證本文所提方法對(duì)衛(wèi)星鐘差預(yù)報(bào)的優(yōu)勢(shì),分別建立了二次多項(xiàng)式模型(QPM)、灰色模型(GM(1,1))、修正指數(shù)曲線法模型(MECM)、自回歸滑動(dòng)平均模型(ARMA)和本文的預(yù)報(bào)模型(New method),對(duì)未來(lái)6 h的衛(wèi)星鐘差進(jìn)行預(yù)報(bào)試驗(yàn),并且將本文預(yù)報(bào)模型預(yù)報(bào)的結(jié)果同其他各模型預(yù)報(bào)的結(jié)果進(jìn)行比較分析.評(píng)價(jià)模型預(yù)報(bào)性能的好壞,采用均方根誤差RMS(計(jì)算公式見(jiàn)式(29))和最大誤差與最小誤差之差的絕對(duì)值,即極差Range(計(jì)算公式見(jiàn)式(30))作為評(píng)價(jià)預(yù)報(bào)結(jié)果的統(tǒng)計(jì)量,去比較分析各模型的預(yù)報(bào)精度及其穩(wěn)定度.其中RMS表征了預(yù)報(bào)結(jié)果的精度,Range表征了算法的穩(wěn)定度.
(29)
(30)
圖5為5種模型的預(yù)報(bào)誤差變化圖,圖6為5種模型的平均預(yù)報(bào)精度RMS和算法穩(wěn)定度Range的統(tǒng)計(jì)圖,表2為5種模型的預(yù)報(bào)誤差的統(tǒng)計(jì)結(jié)果,由圖5~6和表2可知:
圖5 6 h預(yù)報(bào)誤差對(duì)比
表2 衛(wèi)星鐘差預(yù)報(bào)誤差統(tǒng)計(jì)
1)在衛(wèi)星鐘差6 h的預(yù)報(bào)中,本文提出的新方法預(yù)報(bào)精度最高.無(wú)論是對(duì)于早期發(fā)射的IIR型衛(wèi)星,還是最近幾年發(fā)射的IIF型衛(wèi)星,均可以達(dá)到亞納秒級(jí)的預(yù)報(bào)精度,可以實(shí)現(xiàn)厘米級(jí)的定位.算法的穩(wěn)定度也均小于其他預(yù)報(bào)模型,具有較高的預(yù)報(bào)穩(wěn)定性.
2)在衛(wèi)星鐘差6 h的預(yù)報(bào)中,MECM預(yù)報(bào)模型的預(yù)報(bào)性能最差,平均預(yù)報(bào)精度達(dá)到了2.18 ns,平均算法穩(wěn)定度達(dá)到了3.88 ns,這進(jìn)一步說(shuō)明了MECM模型比較適合于類(lèi)指數(shù)變化趨勢(shì)的衛(wèi)星鐘差預(yù)報(bào),而不適合大致呈單調(diào)遞增或者單調(diào)遞減的衛(wèi)星鐘差的預(yù)報(bào);其次是QPM模型,平均預(yù)報(bào)精度為2.01 ns,平均算法穩(wěn)定度為2.36 ns;而GM(1,1)模型和ARMA模型的預(yù)報(bào)性能大致相當(dāng),其平均預(yù)報(bào)精度分別為0.75 ns和0.62 ns,平均算法穩(wěn)定度分別為1.24 ns和1.18 ns.
3)粒子群優(yōu)化加權(quán)灰色回歸組合的自適應(yīng)衛(wèi)星鐘差預(yù)報(bào)模型的精度明顯優(yōu)于其他4種模型,6 h的平均預(yù)報(bào)精度為0.42 ns,平均算法穩(wěn)定度為0.87 ns,相對(duì)于目前正在衛(wèi)星上使用的二次多項(xiàng)式模型的平均預(yù)報(bào)精度提高了79.10%,平均算法穩(wěn)定度提高了63.10%.
4)采用本文的預(yù)報(bào)方法,四顆星載鐘的預(yù)報(bào)精度均達(dá)到了亞納秒量級(jí)的預(yù)報(bào)精度.由于星載鐘受發(fā)射時(shí)間的長(zhǎng)短、設(shè)備老化、鐘差變化趨勢(shì)以及所處環(huán)境等復(fù)雜因素的影響,所以預(yù)報(bào)誤差存在一定程度上的差異.
本文提出的粒子群算法優(yōu)化指數(shù)函數(shù)和線性函數(shù)逼近加權(quán)灰色回歸組合的自適應(yīng)衛(wèi)星鐘差預(yù)報(bào)模型,在建模之前先對(duì)鐘差的質(zhì)量進(jìn)行檢測(cè),剔除異常鐘差數(shù)據(jù)用分段線性插值法將缺失的鐘差數(shù)據(jù)補(bǔ)齊,然后進(jìn)行平滑性檢測(cè),對(duì)于非平滑序列進(jìn)行平滑處理,之后用指數(shù)函數(shù)和線性函數(shù)去逼近灰色模型的一次累加序列后,采用變權(quán)加權(quán)法進(jìn)行組合預(yù)報(bào).針對(duì)精度遞增因子難以確定的問(wèn)題,引入了粒子群優(yōu)化算法進(jìn)行自適應(yīng)尋優(yōu)最佳的精度遞增因子,最后進(jìn)行了6 h的鐘差預(yù)報(bào).在6 h的鐘差預(yù)報(bào)中,所提方法的平均預(yù)報(bào)精度和平均算法穩(wěn)定度均可以達(dá)到亞納秒量級(jí),能夠滿(mǎn)足精密單點(diǎn)定位對(duì)衛(wèi)星鐘差精度的要求,為衛(wèi)星鐘差預(yù)報(bào)提供了一種新的方案,但是該方法仍然沒(méi)有克服誤差累積的現(xiàn)象,所以還有待進(jìn)一步研究如何實(shí)現(xiàn)衛(wèi)星鐘差高精度的預(yù)報(bào).