楊承午，劉志平，徐永明

(1.中國(guó)礦業(yè)大學(xué) 自然資源部國(guó)土環(huán)境與災(zāi)害監(jiān)測(cè)重點(diǎn)實(shí)驗(yàn)室，江蘇 徐州 221116；2.中國(guó)礦業(yè)大學(xué) 環(huán)境與測(cè)繪學(xué)院，江蘇 徐州 221116)

0 引 言

星載原子鐘鐘差的精確測(cè)定是提升GNSS系統(tǒng)定位、導(dǎo)航和授時(shí)功能亟待解決的問題之一。然而，IGS服務(wù)組織及各數(shù)據(jù)分析中心提供的精密鐘差數(shù)據(jù)均存在時(shí)延滯后性，不能滿足用戶實(shí)時(shí)性的需求[1]。為此，國(guó)內(nèi)外學(xué)者在衛(wèi)星鐘差預(yù)報(bào)方面開展了廣泛研究，常見的鐘差預(yù)報(bào)模型有二次多項(xiàng)式模型、灰色系統(tǒng)模型[2-4]、最小二乘支持向量機(jī)模型[5]、Kalman濾波模型[6-8]、神經(jīng)網(wǎng)絡(luò)模型[9-10]和相應(yīng)組合模型[11]。文獻(xiàn)[2]提出將鐘差變化看作一個(gè)灰色系統(tǒng)，將灰色模型引入鐘差預(yù)報(bào)領(lǐng)域，提高了鐘差預(yù)報(bào)精度。現(xiàn)有文獻(xiàn)對(duì)灰色模型的改進(jìn)大多從以下方面入手：一是對(duì)模型參數(shù)中的發(fā)展系數(shù)a進(jìn)行優(yōu)化;二是對(duì)模型參數(shù)中的灰色作用量u進(jìn)行優(yōu)化;三是對(duì)時(shí)間響應(yīng)函數(shù)初始值的優(yōu)化選定[12-13]。這些研究從不同角度對(duì)灰色模型進(jìn)行優(yōu)化，并在實(shí)際應(yīng)用中取得了良好的效果，但這些優(yōu)化都是針對(duì)單個(gè)參數(shù)，并沒有進(jìn)行整體參數(shù)尋優(yōu)。鑒于此，為了獲得更好的建模和預(yù)測(cè)效果，本文引入Levenberg-Marquardt算法[14](以下簡(jiǎn)稱L-M算法)，對(duì)灰色模型參數(shù)進(jìn)行迭代優(yōu)化的改進(jìn)方法。該改進(jìn)方法將灰色模型的模型參數(shù)作為初始值，通過L-M算法對(duì)參數(shù)進(jìn)行迭代尋優(yōu)，尋找出與原始數(shù)據(jù)最為擬合的發(fā)展系數(shù)、灰色作用量和時(shí)間響應(yīng)函數(shù)初始值。

1 鐘差預(yù)報(bào)模型

1.1 灰色GM(1，1)模型

星載原子鐘極易受到外界及其自身因素的影響，穩(wěn)定性差異較大，從而很難掌握其復(fù)雜細(xì)致的變化規(guī)律，鐘差數(shù)據(jù)符合灰色系統(tǒng)理論的建模預(yù)測(cè)特點(diǎn)[2]。

1.1.1 原始鐘差數(shù)據(jù)序列GM(1，1)

若X(0)為非負(fù)準(zhǔn)光滑序列，則其一階累加生成序列X(1)具有準(zhǔn)指數(shù)規(guī)律。一階累加生成序列的背景值Z(1)(k)=λx(1)(k)+(1-λ)x(1)·(k-1)，背景值生成因子λ∈[0，1](一般取為0.5)。則可以構(gòu)造離散化的灰色差分方程

(1)

(2)

式中：k為參與建模的原始數(shù)據(jù)序列長(zhǎng)度；p為預(yù)報(bào)的數(shù)據(jù)序列長(zhǎng)度。

1.1.2 歷元間差分鐘差數(shù)據(jù)序列GM(1，1)

若原始數(shù)據(jù)序列本身呈準(zhǔn)指數(shù)規(guī)律，取原始數(shù)據(jù)序列為X(1)，一階累減生成序列(1-IAGO)為X(0)，構(gòu)建灰色GM(1，1)模型。這種模型構(gòu)建方法可以在一定程度上消除原始鐘差序列趨勢(shì)項(xiàng)的影響，得到一組有效數(shù)字位數(shù)減少、便于進(jìn)行預(yù)處理的數(shù)據(jù)序列，當(dāng)差分?jǐn)?shù)據(jù)序列波動(dòng)較為平穩(wěn)時(shí)，灰色GM(1，1)模型具有更高的預(yù)報(bào)精度[15]。

1.2 基于L-M算法的GM(1，1)模型參數(shù)優(yōu)化

(4)

(7)

L-M算法是Gauss-Newton法與梯度下降法的結(jié)合，保證了全局最優(yōu)收斂，其算法原理如下。

從病犬主人口述，總結(jié)出造成犬患細(xì)小病毒病的幾個(gè)細(xì)節(jié)問題：提早出去遛狗，未打完所有免疫之前不要外出；喂大量肉食；免疫失敗，病痊愈后15～20 d后重新免疫；將食物扔在地上喂食；遛狗時(shí)沒有好的習(xí)慣，小狗亂聞、挖坑吃土之類。綜上所述，關(guān)于犬細(xì)小病毒的診治預(yù)防還需要各方面的積極配合，以降低本地區(qū)該病的發(fā)病率及死亡率。

將式(4)改為函數(shù)表達(dá)式f(t，α0，α1，β，x)=α0eα1(t-1)+β-x，其中t為數(shù)據(jù)序列時(shí)間節(jié)點(diǎn)，x為建模數(shù)據(jù)序列值，根據(jù)泰勒定理可得其一階展開式

f(t，α0+h(1)，α1+h(2)，β+h(3)，x)?

(8)

進(jìn)一步，可建立f(t，α0，α1，β，x)的均值方差準(zhǔn)則函數(shù)L(h)，即

L(h)=(h)T(h)/2=fTf+hTJTf/2+

hTJTJh/2,

(9)

h=-(JTJ+μI)-1JTf,

(10)

式中，μ為阻尼系數(shù)，協(xié)方差矩陣(JTJ+μI)對(duì)稱正定，確保h代表f(x)的下降方向。

式(10)中，阻尼系數(shù)μ=τmax {diag(A)}，A=JTJ，為設(shè)計(jì)矩陣，τ為阻尼系數(shù)變換因子。在迭代開始時(shí)，τ由用戶根據(jù)實(shí)際最優(yōu)化問題設(shè)定一個(gè)較小正值；在迭代過程中通過增益比ρ[15]進(jìn)行修正。

ρ=(2hTJTf+hTJTJh)hTμh-hTJTf，

(11)

τ(k+1)=

分析式(11)，(12)可知，由于hTh和-hTJTf=hT(JTJ+μI)h均為非負(fù)，則可以保證增益比ρ的分母一定非負(fù)。若ρ為正值，此時(shí)應(yīng)減小阻尼系數(shù)μ以增加迭代步長(zhǎng)，一般取a=1/3，使下次迭代更接近Gauss-newton算法，發(fā)揮其局部收斂特性且加速收斂速度；若ρ為負(fù)值，此時(shí)應(yīng)增大阻尼系數(shù)μ,以減小迭代步長(zhǎng)，一般取b=2，使下次迭代更接近梯度下降算法，發(fā)揮其全局特性并保證全局收斂[15]。施加阻尼系數(shù)μ既可以控制下降方向又可以控制下降步長(zhǎng)，保證迭代的全局最優(yōu)收斂和較高的收斂速度。

基于L-M算法的基本原理與公式闡述，以下對(duì)其在GM(1，1)模型鐘差預(yù)報(bào)應(yīng)用中的流程進(jìn)行簡(jiǎn)單描述，步驟如下。

(3)參數(shù)回代：將優(yōu)化后的參數(shù)回代時(shí)間響應(yīng)函數(shù)，求解鐘差序列的擬合與預(yù)測(cè)值。

2 GPS星載原子鐘LM-GM預(yù)報(bào)實(shí)驗(yàn)

為驗(yàn)證本文提出的基于L-M算法對(duì)GM(1，1)參數(shù)估計(jì)改進(jìn)的有效性，采用IGS發(fā)布的2018年5月6—7日采樣時(shí)間間隔15 min的精密衛(wèi)星鐘差數(shù)據(jù)進(jìn)行仿真分析。該時(shí)間段內(nèi)在軌運(yùn)行的星載原子鐘共分為4種類型，如表1所示，隨機(jī)選取每類衛(wèi)星中的一顆，本文選取PRN06、PRN07、PRN22和PRN24的衛(wèi)星鐘差數(shù)據(jù)進(jìn)行實(shí)驗(yàn)。

表1 不同原子鐘類型及對(duì)應(yīng)的衛(wèi)星編號(hào)Tab.1 Types of different atomic clock and corresponding satellite number

IGS服務(wù)器提供的超快速星歷產(chǎn)品更新周期為6 h，快速星歷產(chǎn)品更新周期為24 h。鑒于此，本文設(shè)計(jì)2種預(yù)報(bào)方案，使用IGS服務(wù)器提供的1 d精密星歷產(chǎn)品構(gòu)建模型。方案一：預(yù)報(bào)6 h的短期鐘差數(shù)據(jù)；方案二：預(yù)報(bào)24 h的中期鐘差數(shù)據(jù)。比較分析2種方案下，本文所提LM-GM優(yōu)化模型與二次多項(xiàng)式模型(QP)和灰色模型(GM)鐘差數(shù)據(jù)預(yù)報(bào)精度的差異，預(yù)報(bào)誤差統(tǒng)計(jì)比較結(jié)果如圖1和圖2所示。

由圖1和圖2可知，隨著預(yù)報(bào)時(shí)間的增加，3種預(yù)報(bào)模型的預(yù)報(bào)誤差均呈現(xiàn)增長(zhǎng)趨勢(shì)。LM-GM優(yōu)化模型相較其余2種模型可以取得更高的預(yù)報(bào)精度，對(duì)于穩(wěn)定性較差的CS原子鐘的鐘差預(yù)報(bào)，GM模型和LM-GM模型可以取得比QP模型更高的預(yù)報(bào)精度。

進(jìn)一步，結(jié)合圖1、圖2和表 2進(jìn)行結(jié)果分析，可以得出如下結(jié)論：

(1)4顆衛(wèi)星原子鐘鐘差數(shù)據(jù)采用3種預(yù)報(bào)模型和2種預(yù)報(bào)方案，均可以滿足理想的預(yù)報(bào)精度要求。本文所提的參數(shù)估計(jì)優(yōu)化方法的灰色模型(LM-GM)在不同方案下的預(yù)報(bào)精度均優(yōu)于二次多項(xiàng)式模型(QP)和灰色模型(GM)的預(yù)報(bào)精度。對(duì)于6 h短期鐘差預(yù)報(bào)方案，LM-GM模型Rb原子鐘的預(yù)報(bào)精度均處于亞納秒量級(jí)，相較二次多項(xiàng)式模型，預(yù)報(bào)精度提高了0.4%～68.7%，相較灰色模型，預(yù)報(bào)精度提高了43.3%～83.3%；LM-GM模型Cs原子鐘的預(yù)報(bào)精度為1.876 6 ns，相較于QP模型的4.745 1 ns和GM模型的2.083 5 ns，預(yù)報(bào)精度提高了60.5%和10%；預(yù)報(bào)精度平均提高了48.5%。對(duì)于24 h中期鐘差預(yù)報(bào)方案，3種模型的預(yù)報(bào)精度均顯著降低，LM-GM模型Rb原子鐘的預(yù)報(bào)精度處于亞納秒至納秒量級(jí)之間，相較于二次多項(xiàng)式模型，預(yù)報(bào)精度提高了3%～61.3%，相較灰色模型，預(yù)報(bào)精度提高了35.2%～82.6%；LM-GM模型Cs原子鐘的預(yù)報(bào)精度為4.871 6 ns，相較于QP模型的17.202 2 ns和GM模型的5.119 4 ns，預(yù)報(bào)精度提高了71.7%和4.9%；預(yù)報(bào)精度平均提高了46.3%。

(2)由于Cs原子鐘短期穩(wěn)定性較差，星載Rb原子鐘的預(yù)報(bào)精度要明顯高于星載Cs原子鐘的預(yù)報(bào)精度。在進(jìn)行鐘差數(shù)據(jù)建模預(yù)測(cè)時(shí)，使用Rb原子鐘鐘差數(shù)據(jù)可以取得更高的預(yù)報(bào)精度。目前GPS衛(wèi)星系統(tǒng)搭載的原子鐘逐漸替換為Rb原子鐘，北斗導(dǎo)航系統(tǒng)全部搭載Rb原子鐘，這樣可以獲取到更為穩(wěn)定的原始鐘差數(shù)據(jù)和更精確的預(yù)報(bào)數(shù)據(jù)，從而在精密單點(diǎn)定位、導(dǎo)航、授時(shí)等應(yīng)用中滿足用戶更高的需求。

圖1 方案一預(yù)報(bào)誤差結(jié)果對(duì)比Fig.1 Prediction error result comparison of scheme one

圖2 方案二預(yù)報(bào)誤差結(jié)果對(duì)比Fig.2 Prediction error result comparison of scheme two

圖3 L-M算法參數(shù)隨迭代過程的變化Fig.3 Parameter variation of L-M algorithm in iteration process

3 結(jié) 語

本文提出了一種基于L-M算法對(duì)灰色模型參數(shù)迭代優(yōu)化方法，該方法將灰色模型參數(shù)估值作為迭代初始值，通過L-M算法尋找出與原始數(shù)據(jù)最優(yōu)擬合的灰色模型參數(shù)估值，利用優(yōu)化的模型參數(shù)求解時(shí)間響應(yīng)函數(shù)，提高了灰色模型的預(yù)報(bào)精度。結(jié)合IGS服務(wù)中心提供的精密鐘差數(shù)據(jù)產(chǎn)品進(jìn)行預(yù)報(bào)實(shí)驗(yàn)，結(jié)果表明，相較傳統(tǒng)預(yù)報(bào)模型，改進(jìn)模型的預(yù)報(bào)精度均得到提升。6 h預(yù)報(bào)方案中，Rb原子鐘的預(yù)報(bào)精度處于亞納秒量級(jí)，Cs原子鐘的預(yù)報(bào)精度優(yōu)于2 ns；24 h預(yù)報(bào)方案中，Rb原子鐘的預(yù)報(bào)精度處于亞納秒至納秒量級(jí)，Cs原子鐘的預(yù)報(bào)精度優(yōu)于5 ns，能夠滿足目前導(dǎo)航定位對(duì)衛(wèi)星鐘差精度的要求。