粒子群優化加權灰色回歸組合的衛星鐘差預報

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

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

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

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

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

1 衛星鐘差數據預處理

1.1 衛星鐘差數據異常探測

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

(1)

式中Median(yi)為衛星鐘差頻率數據yi的中間數.

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

1.2 衛星鐘差數據平滑性檢測

由于星載原子鐘受自身復雜的時頻特性和外界復雜環境的影響,穩定性存在較大差異的特點,而衛星鐘差是否平滑對于模型的預報性能有很大影響,所以對觀測的鐘差進行平滑性檢測就顯得很有必要.

若ρ(i)滿足:

(2)

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

1.3 三點平滑處理

由于衛星鐘差數據存在系統噪聲,采用三點平滑法處理,可以提高衛星鐘差數據序列的平滑度.三點平滑法是一種取算術平均的方法,它通過重新分配待處理的衛星鐘差與前后鐘差數據的權值,以加強待處理鐘差數據的權重,從而減小衛星鐘差的波動性,增強待處理的衛星鐘差與前后鐘差數據間的聯系.其過程如下:

設有一組衛星鐘差數據序列為

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

(3)

平滑處理位于衛星鐘差中間的數據為

(4)

而對于鐘差數據序列兩端的數據平滑處理的計算為

(5)

經式(4)和(5)處理后,得到一組噪聲減小的鐘差數據序列,然后建模預報.

2 基于PSO優化的灰色加權回歸組合的自適應預報模型

2.1 加權灰色回歸組合模型

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

設X(0)={x(0)(1),x(0)(2),…,x(0)(n)}為不同時刻的衛星鐘差,經過中位數法探測剔除異常鐘差數據,然后采用分段線性插值法把缺失的鐘差數據補齊,得到一組新的鐘差數據序列為X(1)={x(1)(1),x(1)(2),…,x(1)(n)}.再使用三點平滑法處理,以得到一組噪聲減小的鐘差數據序列為X(2)={x(2)(1),x(2)(2),…,x(2)(n)},然后對X(2)這組鐘差數據序列建立灰色模型,則這組數據的一次累加鐘差數據序列為X(3)={x(3)(1),x(3)(2),…,x(3)(n)}的灰色預報模型可表示為

(6)

記參數為A=[au]T,根據最小二乘法可得

(7)

式中

將式(7)代入式(6)即可預報出鐘差的一次累加序列，然后還原即可得到原始鐘差的預報值.

分析可以看出式(6)為以下形式

(8)

由于上式為指數函數和線性函數的疊加形式,下面采用指數曲線y=aex與線性曲線y=αx+β這兩種模型的和來逼近鐘差序列X(2)的一次累加生成的序列X(3)的預報值,如下

(9)

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

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

(10)

設:

(11)

則:

(12)

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

(13)

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

(14)

當m=1時:

(15)

當m=2時:

(16)

直到當m=n-3時:

(17)

(18)

(19)

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

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

(20)

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

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

(21)

式中

精度遞增因子R的選取是該算法預報性能好壞的關鍵,本文采用粒子群優化算法(Particle Swarm Optimization,PSO)對精度遞增因子R進行自適應尋優.

2.2 自適應精度遞增因子的求解算法

PSO算法是由美國兩位科學家受鳥類群體尋找食物的行為得到啟發而提出的一種智能優化算法.因為PSO算法收斂速度比較快,調整參數簡單,所以得到了廣泛的應用.因此,本文選用PSO算法對上述精度遞增因子R進行自適應尋優,PSO算法原理[15]如下:

假設在一個目標搜索空間,其維數為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個粒子迄今為止搜索到的最優位置稱為個體極值,記為

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

(24)

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

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

(25)

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

(26)

式中:c1和c2為加速常數,一般取c1=c2=2;r1和r2為[0,1]范圍內的均勻隨機數;wi為慣性權重,wmax和wmin為預先確定的最大和最小慣性權重,一般wmax取0.9,wmin取0.1;Gmax為最大迭代次數;Gi為當前迭代次數.

采用PSO算法對精度遞增因子R進行自適應尋優的具體步驟如下:

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

(27)

式中:yi為預報的鐘差,yip為實際的觀測鐘差.

2)計算各粒子的適應度函數值fi,然后再比較它們的適應度函數值fi的大小.如果粒子的當前適應度值fi小于該粒子的歷史最優適應度值fibest,則把fibest賦給fi;同理若fibest小于當前粒子種群的最優適應度值fgbest,則把fgbest賦給fibest.

3)根據式(26)更新粒子的速度、位置和慣性權重,然后判斷是否達到了最大迭代次數Gmax.若達到,則迭代結束,所得到的最優位置即為精度遞增因子R的值;若未達到,則繼續迭代直到迭代結束.此算法的具體流程見圖1.

圖1 粒子群優化加權灰色回歸組合自適應預報算法流程

3 試驗與分析

3.1 試驗數據來源

為了驗證本文預報方法的有效性和可行性,采用IGS服務器(ftp://igs.ign.fr/pub/igs/)上發布的GPS 2 000周,間隔時長為15 min的SP3格式的事后精密衛星鐘差數據進行預報試驗.在該時間段內在軌的GPS衛星有31顆,其星載原子鐘有BLOCK IIR-Rb鐘、BLOCK IIR-M-Rb鐘、BLOCK IIF-Rb鐘和BLOCK IIF-Cs四種類型.由于北斗衛星導航系統的星載鐘與GPS基本一致,且北斗二代衛星導航系統均搭載的為銣原子鐘,為使研究結果能為我國北斗衛星導航系統在鐘差預報方面提供一些參考,所以選取了GPS IIF-Rb PRN03、GPS IIR-M-Rb PRN12、GPS IIR-Rb PRN14和GPS IIR-M-Rb PRN17四顆搭載銣原子鐘的衛星的鐘差數據進行預報試驗.截止2019年10月16日,它們的相關信息見表1.

表1 選擇的衛星相關信息

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

圖2 原始觀測衛星鐘差變化趨勢

由于衛星鐘差的有效位數比較多，使得原始觀測鐘差數據異常點容易被掩蓋,而異常鐘差在其對應的頻率數據中表現為顯著的峰值點,所以先將衛星鐘差轉換為相應的頻率數據后更容易對異常鐘差進行探測,具體轉換公式如下

(28)

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

PRN03、PRN12、PRN14和PRN17號衛星相應的頻率數據變化見圖3,經中位數法探測發現,在此時間段內,PRN03、PRN14和PRN17三顆衛星存在異常的鐘差數據,將探測的異常鐘差數據剔除后采用分段線性插值法將缺失的數據補齊,補齊后的變化見圖4,然后還原為鐘差的相位數據進行建模預報.

3.2 預報結果與分析

IGS發布的超快速衛星星歷產品的更新時間約為6 h,鑒于此,為了驗證本文所提方法對衛星鐘差預報的優勢,分別建立了二次多項式模型(QPM)、灰色模型(GM(1,1))、修正指數曲線法模型(MECM)、自回歸滑動平均模型(ARMA)和本文的預報模型(New method),對未來6 h的衛星鐘差進行預報試驗,并且將本文預報模型預報的結果同其他各模型預報的結果進行比較分析.評價模型預報性能的好壞,采用均方根誤差RMS(計算公式見式(29))和最大誤差與最小誤差之差的絕對值,即極差Range(計算公式見式(30))作為評價預報結果的統計量,去比較分析各模型的預報精度及其穩定度.其中RMS表征了預報結果的精度,Range表征了算法的穩定度.

(29)

(30)

圖5為5種模型的預報誤差變化圖,圖6為5種模型的平均預報精度RMS和算法穩定度Range的統計圖，表2為5種模型的預報誤差的統計結果,由圖5～6和表2可知:

圖5 6 h預報誤差對比

表2 衛星鐘差預報誤差統計

1)在衛星鐘差6 h的預報中,本文提出的新方法預報精度最高.無論是對于早期發射的IIR型衛星,還是最近幾年發射的IIF型衛星,均可以達到亞納秒級的預報精度,可以實現厘米級的定位.算法的穩定度也均小于其他預報模型,具有較高的預報穩定性.

2)在衛星鐘差6 h的預報中,MECM預報模型的預報性能最差,平均預報精度達到了2.18 ns,平均算法穩定度達到了3.88 ns,這進一步說明了MECM模型比較適合于類指數變化趨勢的衛星鐘差預報,而不適合大致呈單調遞增或者單調遞減的衛星鐘差的預報;其次是QPM模型,平均預報精度為2.01 ns,平均算法穩定度為2.36 ns;而GM(1,1)模型和ARMA模型的預報性能大致相當,其平均預報精度分別為0.75 ns和0.62 ns,平均算法穩定度分別為1.24 ns和1.18 ns.

3)粒子群優化加權灰色回歸組合的自適應衛星鐘差預報模型的精度明顯優于其他4種模型,6 h的平均預報精度為0.42 ns,平均算法穩定度為0.87 ns,相對于目前正在衛星上使用的二次多項式模型的平均預報精度提高了79.10%,平均算法穩定度提高了63.10%.

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

4 結論

本文提出的粒子群算法優化指數函數和線性函數逼近加權灰色回歸組合的自適應衛星鐘差預報模型,在建模之前先對鐘差的質量進行檢測,剔除異常鐘差數據用分段線性插值法將缺失的鐘差數據補齊,然后進行平滑性檢測,對于非平滑序列進行平滑處理,之后用指數函數和線性函數去逼近灰色模型的一次累加序列后,采用變權加權法進行組合預報.針對精度遞增因子難以確定的問題,引入了粒子群優化算法進行自適應尋優最佳的精度遞增因子,最后進行了6 h的鐘差預報.在6 h的鐘差預報中,所提方法的平均預報精度和平均算法穩定度均可以達到亞納秒量級,能夠滿足精密單點定位對衛星鐘差精度的要求,為衛星鐘差預報提供了一種新的方案,但是該方法仍然沒有克服誤差累積的現象,所以還有待進一步研究如何實現衛星鐘差高精度的預報.