于 燁，黃 默，段 濤，王長元，胡 銳

(1.中國科學(xué)院 微電子研究所，北京 100029； 2.中國科學(xué)院大學(xué)，北京 100049；3.中國科學(xué)院大學(xué) 微電子學(xué)院，北京 100049)

衛(wèi)星鐘差預(yù)報(bào)是導(dǎo)航定位中一項(xiàng)至關(guān)重要的工作之一[1-2].星載原子鐘因其自身的性能以及所處空間復(fù)雜環(huán)境的影響,鐘源偏差需要與地面站的原子鐘進(jìn)行比對并校準(zhǔn)以保持比較精確的時間信息,但是出于某些原因星載原子鐘無法與地面站保持比對的時候就需要采用已有的鐘差信息進(jìn)行預(yù)報(bào).在這段時間內(nèi),衛(wèi)星以預(yù)報(bào)的星歷信息進(jìn)行自主完成軌道確定和廣播星歷的播發(fā),這對于衛(wèi)星的自主生存能力具有重要影響[3-4],所以如何提高衛(wèi)星鐘差預(yù)報(bào)的精度和穩(wěn)定度一直以來都是國內(nèi)外研究的熱點(diǎn)問題之一.

長久以來,在衛(wèi)星鐘差預(yù)報(bào)方面,國內(nèi)外許多學(xué)者進(jìn)行了多角度、多方位和深層次的研究,并取得了一系列重要研究成果[5-7].其中,灰色系統(tǒng)預(yù)報(bào)模型只需要少量的鐘差數(shù)據(jù)建模就可以預(yù)報(bào),并且需要數(shù)據(jù)平滑且呈類指數(shù)變化,對于數(shù)據(jù)呈非指數(shù)變化的鐘差預(yù)報(bào)效果較差等特點(diǎn).文獻(xiàn)[8]提出采用最小一乘法準(zhǔn)則去估計(jì)灰色模型的參數(shù),以克服在鐘差波動較大的情況下最小二乘法準(zhǔn)則估計(jì)參數(shù)預(yù)報(bào)精度的不足,并且通過預(yù)報(bào)試驗(yàn)證明了算法的有效性.文獻(xiàn)[9]提出了一種先用對數(shù)平滑處理的方法對原始鐘差數(shù)據(jù)序列進(jìn)行平滑處理,再使用自適應(yīng)雙子群改進(jìn)粒子群算法來優(yōu)化灰色模型的發(fā)展系數(shù)和灰作用量的預(yù)報(bào)策略,從而增強(qiáng)了灰色模型的泛化能力.文獻(xiàn)[10]提出先對衛(wèi)星鐘差第一天偶數(shù)整時的值采用灰色模型進(jìn)行建模,然后利用已建立的模型預(yù)報(bào)第二天相應(yīng)時間的鐘差,對預(yù)報(bào)得到的鐘差進(jìn)行拉格朗日插值建模,通過內(nèi)插得到其他時刻的預(yù)報(bào)值,以此來提高模型的預(yù)報(bào)能力.文獻(xiàn)[11]分析了灰色模型發(fā)展系數(shù)和灰作用量對預(yù)報(bào)精度的影響,之后通過微調(diào)這兩個參數(shù)來提高灰色模型的預(yù)報(bào)性能,但是針對鐘差數(shù)據(jù)序列呈現(xiàn)出不同變化趨勢時,沒有給出有效的參數(shù)優(yōu)化方法.

從以上的分析可知,研究者們主要是在優(yōu)化灰色模型的兩個參數(shù)方面對模型進(jìn)行改進(jìn),沒有從灰色模型的一次累加生成序列的預(yù)報(bào)模型方面進(jìn)行其他處理.基于此,本文提出了一種基于粒子群算法優(yōu)化指數(shù)函數(shù)和線性函數(shù)逼近加權(quán)灰色回歸組合的自適應(yīng)衛(wèi)星鐘差預(yù)報(bào)方法.在建模之前,考慮到衛(wèi)星鐘差鐘跳頻繁的現(xiàn)象,首先采用中位數(shù)法對衛(wèi)星鐘差進(jìn)行異常鐘差探測,將異常鐘差剔除后采用分段線性插值法將缺失的鐘差補(bǔ)齊.然后考慮到衛(wèi)星鐘差存在系統(tǒng)噪聲,采用三點(diǎn)平滑法對鐘差數(shù)據(jù)進(jìn)行平滑處理后建立了粒子群算法優(yōu)化指數(shù)函數(shù)和線性函數(shù)逼近加權(quán)灰色回歸組合的自適應(yīng)衛(wèi)星鐘差預(yù)報(bào)模型,并結(jié)合4種典型變化趨勢的衛(wèi)星鐘差序列,采用IGS服務(wù)器上發(fā)布的事后精密衛(wèi)星鐘差產(chǎn)品作為試驗(yàn)數(shù)據(jù),進(jìn)行了預(yù)報(bào)試驗(yàn).

1 衛(wèi)星鐘差數(shù)據(jù)預(yù)處理

1.1 衛(wèi)星鐘差數(shù)據(jù)異常探測

由于衛(wèi)星鐘差數(shù)據(jù)跳變頻繁,跳變的鐘差數(shù)據(jù)很不利于建模,所以在建模之前對衛(wèi)星鐘差數(shù)據(jù)的質(zhì)量進(jìn)行檢測就顯得很重要.本文采用時效性和抗差性較好的中位數(shù)法(Median Absolute Deviation, MAD)進(jìn)行異常鐘差探測[12],其表達(dá)式為

(1)

式中Median(yi)為衛(wèi)星鐘差頻率數(shù)據(jù)yi的中間數(shù).

當(dāng)衛(wèi)星鐘差的頻率yi>(m+n·MAD)或yi<(m+n·MAD)時,可以判斷該點(diǎn)為異常鐘差.將探測到的異常鐘差剔除后,采用分段線性插值法把缺失的鐘差數(shù)據(jù)補(bǔ)齊.

1.2 衛(wèi)星鐘差數(shù)據(jù)平滑性檢測

由于星載原子鐘受自身復(fù)雜的時頻特性和外界復(fù)雜環(huán)境的影響,穩(wěn)定性存在較大差異的特點(diǎn),而衛(wèi)星鐘差是否平滑對于模型的預(yù)報(bào)性能有很大影響,所以對觀測的鐘差進(jìn)行平滑性檢測就顯得很有必要.

若ρ(i)滿足:

(2)

則稱鐘差數(shù)據(jù)序列X(1)為平滑序列;反之,稱鐘差數(shù)據(jù)序列X(1)為非平滑序列.

1.3 三點(diǎn)平滑處理

由于衛(wèi)星鐘差數(shù)據(jù)存在系統(tǒng)噪聲,采用三點(diǎn)平滑法處理,可以提高衛(wèi)星鐘差數(shù)據(jù)序列的平滑度.三點(diǎn)平滑法是一種取算術(shù)平均的方法,它通過重新分配待處理的衛(wèi)星鐘差與前后鐘差數(shù)據(jù)的權(quán)值,以加強(qiáng)待處理鐘差數(shù)據(jù)的權(quán)重,從而減小衛(wèi)星鐘差的波動性,增強(qiáng)待處理的衛(wèi)星鐘差與前后鐘差數(shù)據(jù)間的聯(lián)系.其過程如下:

設(shè)有一組衛(wèi)星鐘差數(shù)據(jù)序列為

X(1)={x(1)(1),x(1)(2),…,x(1)(n)}，

(3)

平滑處理位于衛(wèi)星鐘差中間的數(shù)據(jù)為

(4)

而對于鐘差數(shù)據(jù)序列兩端的數(shù)據(jù)平滑處理的計(jì)算為

(5)

經(jīng)式(4)和(5)處理后,得到一組噪聲減小的鐘差數(shù)據(jù)序列,然后建模預(yù)報(bào).

2 基于PSO優(yōu)化的灰色加權(quán)回歸組合的自適應(yīng)預(yù)報(bào)模型

2.1 加權(quán)灰色回歸組合模型

灰色模型是指信息不完全知道的系統(tǒng),具體是指它的部分信息是已知的、其余的信息為未知的一種系統(tǒng).其最常見的灰色系統(tǒng)是GM(1,1)模型,它由一個包含單變量的一階微分方程所構(gòu)成,可以實(shí)現(xiàn)對自身數(shù)據(jù)的預(yù)測[13].

設(shè)X(0)={x(0)(1),x(0)(2),…,x(0)(n)}為不同時刻的衛(wèi)星鐘差,經(jīng)過中位數(shù)法探測剔除異常鐘差數(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)},然后對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+β這兩種模型的和來逼近鐘差序列X(2)的一次累加生成的序列X(3)的預(yù)報(bào)值,如下

(9)

式中,參數(shù)v和參數(shù)C1、C2、C3為待求的參數(shù).

為確定式(9)中的4個未知參數(shù)，構(gòu)造如下序列:

(10)

設(shè):

(11)

則:

(12)

用式(12)除以式(11)得:

(13)

將式(13)兩邊取對數(shù)得:

(14)

當(dāng)m=1時:

(15)

當(dāng)m=2時:

(16)

直到當(dāng)m=n-3時:

(17)

(18)

(19)

由于離預(yù)報(bào)值近一點(diǎn)的鐘差對預(yù)報(bào)的鐘差影響大一點(diǎn),離預(yù)報(bào)值遠(yuǎn)一點(diǎn)的鐘差對預(yù)報(bào)的鐘差影響小一點(diǎn),下面采用變權(quán)加權(quán)法來解式(19).變權(quán)加權(quán)法相對于普通加權(quán)法最大的不同是,變權(quán)加權(quán)法對鐘差時間序列的可靠性成比例的分別賦予不同的權(quán)值[14],具體加權(quán)方法如下:

Pk=Ri-1,i=1,2,…,n.

(20)

式中:R稱為精度遞增因子且R∈[1,2].

綜合式(19)和式(20),利用最小二乘法可得

(21)

式中

精度遞增因子R的選取是該算法預(yù)報(bào)性能好壞的關(guān)鍵,本文采用粒子群優(yōu)化算法(Particle Swarm Optimization,PSO)對精度遞增因子R進(jìn)行自適應(yīng)尋優(yōu).

2.2 自適應(yīng)精度遞增因子的求解算法

PSO算法是由美國兩位科學(xué)家受鳥類群體尋找食物的行為得到啟發(fā)而提出的一種智能優(yōu)化算法.因?yàn)镻SO算法收斂速度比較快,調(diào)整參數(shù)簡單,所以得到了廣泛的應(yīng)用.因此,本文選用PSO算法對上述精度遞增因子R進(jìn)行自適應(yīng)尋優(yōu),PSO算法原理[15]如下:

假設(shè)在一個目標(biāo)搜索空間,其維數(shù)為D維且有N個粒子組成一個群落,其中第i個粒子可表示為一個D維向量

Xi=(xi1,xi2,…,xiD),i=1,2,…,N.

(22)

第i個粒子的“飛行”速度也是一個D維向量,記為

Vi=(vi1,vi2,…,viD),i=1,2,…,N.

(23)

第i個粒子迄今為止搜索到的最優(yōu)位置稱為個體極值,記為

pbest=(pi1,pi2,…,piD),i=1,2,…,N.

(24)

整個粒子群迄今為止搜索到的最優(yōu)位置為全局最優(yōu)值,記為

gbest=(pg1,pg2,…,pgD).

(25)

在未找到這兩個最優(yōu)值時,粒子根據(jù)如下的公式來更新自己的速度和位置:

(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算法對精度遞增因子R進(jìn)行自適應(yīng)尋優(yōu)的具體步驟如下:

1)初始化粒子群的速度、位置、種群規(guī)模、最大迭代次數(shù)、最大最小慣性權(quán)重、訓(xùn)練樣本數(shù)m和加速常數(shù),同時定義粒子群的適應(yīng)度函數(shù)為

(27)

式中:yi為預(yù)報(bào)的鐘差,yip為實(shí)際的觀測鐘差.

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é)束.此算法的具體流程見圖1.

圖1 粒子群優(yōu)化加權(quán)灰色回歸組合自適應(yīng)預(yù)報(bào)算法流程

3 試驗(yàn)與分析

3.1 試驗(yàn)數(shù)據(jù)來源

為了驗(yàn)證本文預(yù)報(bào)方法的有效性和可行性,采用IGS服務(wù)器(ftp://igs.ign.fr/pub/igs/)上發(fā)布的GPS 2 000周,間隔時長為15 min的SP3格式的事后精密衛(wèi)星鐘差數(shù)據(jù)進(jìn)行預(yù)報(bào)試驗(yàn).在該時間段內(nèi)在軌的GPS衛(wèi)星有31顆,其星載原子鐘有BLOCK IIR-Rb鐘、BLOCK IIR-M-Rb鐘、BLOCK IIF-Rb鐘和BLOCK IIF-Cs四種類型.由于北斗衛(wèi)星導(dǎo)航系統(tǒng)的星載鐘與GPS基本一致,且北斗二代衛(wèi)星導(dǎo)航系統(tǒng)均搭載的為銣原子鐘,為使研究結(jié)果能為我國北斗衛(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)信息見表1.

表1 選擇的衛(wèi)星相關(guān)信息

這四顆衛(wèi)星第2 000周第0 d的事后精密鐘差時間序列的變化情況如圖2所示,其中PRN03號衛(wèi)星的鐘差時間序列呈正值單調(diào)遞增趨勢,PRN12號衛(wèi)星的鐘差時間序列呈正值單調(diào)遞減趨勢,PRN14號衛(wèi)星的鐘差時間序列呈負(fù)值遞減趨勢,PRN17號衛(wèi)星的鐘差時間序列呈負(fù)值單調(diào)遞增趨勢,具有充分的代表性.

圖2 原始觀測衛(wèi)星鐘差變化趨勢

由于衛(wèi)星鐘差的有效位數(shù)比較多，使得原始觀測鐘差數(shù)據(jù)異常點(diǎn)容易被掩蓋,而異常鐘差在其對應(yīng)的頻率數(shù)據(jù)中表現(xiàn)為顯著的峰值點(diǎn),所以先將衛(wèi)星鐘差轉(zhuǎn)換為相應(yīng)的頻率數(shù)據(jù)后更容易對異常鐘差進(jìn)行探測,具體轉(zhuǎn)換公式如下

(28)

式中:yi表示衛(wèi)星鐘差的頻率數(shù)據(jù),xi+1和xi分別表 示第i+1和i歷元的衛(wèi)星鐘差值(相位數(shù)據(jù)),τ0表示采樣間隔.

PRN03、PRN12、PRN14和PRN17號衛(wèi)星相應(yīng)的頻率數(shù)據(jù)變化見圖3,經(jīng)中位數(shù)法探測發(fā)現(xiàn),在此時間段內(nèi),PRN03、PRN14和PRN17三顆衛(wèi)星存在異常的鐘差數(shù)據(jù),將探測的異常鐘差數(shù)據(jù)剔除后采用分段線性插值法將缺失的數(shù)據(jù)補(bǔ)齊,補(bǔ)齊后的變化見圖4,然后還原為鐘差的相位數(shù)據(jù)進(jìn)行建模預(yù)報(bào).

3.2 預(yù)報(bào)結(jié)果與分析

IGS發(fā)布的超快速衛(wèi)星星歷產(chǎn)品的更新時間約為6 h,鑒于此,為了驗(yàn)證本文所提方法對衛(wèi)星鐘差預(yù)報(bào)的優(yōu)勢,分別建立了二次多項(xiàng)式模型(QPM)、灰色模型(GM(1,1))、修正指數(shù)曲線法模型(MECM)、自回歸滑動平均模型(ARMA)和本文的預(yù)報(bào)模型(New method),對未來6 h的衛(wèi)星鐘差進(jìn)行預(yù)報(bào)試驗(yàn),并且將本文預(yù)報(bào)模型預(yù)報(bào)的結(jié)果同其他各模型預(yù)報(bào)的結(jié)果進(jìn)行比較分析.評價模型預(yù)報(bào)性能的好壞,采用均方根誤差RMS(計(jì)算公式見式(29))和最大誤差與最小誤差之差的絕對值,即極差Range(計(jì)算公式見式(30))作為評價預(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)誤差對比

表2 衛(wèi)星鐘差預(yù)報(bào)誤差統(tǒng)計(jì)

1)在衛(wèi)星鐘差6 h的預(yù)報(bào)中,本文提出的新方法預(yù)報(bào)精度最高.無論是對于早期發(fā)射的IIR型衛(wèi)星,還是最近幾年發(fā)射的IIF型衛(wèi)星,均可以達(dá)到亞納秒級的預(yù)報(bào)精度,可以實(shí)現(xiàn)厘米級的定位.算法的穩(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)一步說明了MECM模型比較適合于類指數(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,相對于目前正在衛(wèi)星上使用的二次多項(xiàng)式模型的平均預(yù)報(bào)精度提高了79.10%,平均算法穩(wěn)定度提高了63.10%.

4)采用本文的預(yù)報(bào)方法,四顆星載鐘的預(yù)報(bào)精度均達(dá)到了亞納秒量級的預(yù)報(bào)精度.由于星載鐘受發(fā)射時間的長短、設(shè)備老化、鐘差變化趨勢以及所處環(huán)境等復(fù)雜因素的影響,所以預(yù)報(bào)誤差存在一定程度上的差異.

4 結(jié)論

本文提出的粒子群算法優(yōu)化指數(shù)函數(shù)和線性函數(shù)逼近加權(quán)灰色回歸組合的自適應(yīng)衛(wèi)星鐘差預(yù)報(bào)模型,在建模之前先對鐘差的質(zhì)量進(jìn)行檢測,剔除異常鐘差數(shù)據(jù)用分段線性插值法將缺失的鐘差數(shù)據(jù)補(bǔ)齊,然后進(jìn)行平滑性檢測,對于非平滑序列進(jìn)行平滑處理,之后用指數(shù)函數(shù)和線性函數(shù)去逼近灰色模型的一次累加序列后,采用變權(quán)加權(quán)法進(jìn)行組合預(yù)報(bào).針對精度遞增因子難以確定的問題,引入了粒子群優(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á)到亞納秒量級,能夠滿足精密單點(diǎn)定位對衛(wèi)星鐘差精度的要求,為衛(wèi)星鐘差預(yù)報(bào)提供了一種新的方案,但是該方法仍然沒有克服誤差累積的現(xiàn)象,所以還有待進(jìn)一步研究如何實(shí)現(xiàn)衛(wèi)星鐘差高精度的預(yù)報(bào).