雷騏瑋,張洪波

國防科技大學(xué) 空天科學(xué)學(xué)院,長沙 410073

1 引言

建設(shè)空間站是載人航天發(fā)展的高級階段,高精度的軌道預(yù)報(bào)是空間站在軌安全穩(wěn)定運(yùn)行和制定合理的軌道控制策略的前提條件?？臻g站一般運(yùn)行在低地球軌道,且面質(zhì)比較大,大氣阻力攝動是空間站軌道預(yù)報(bào)的主要誤差源[1]。

計(jì)算低軌航天器所受的大氣阻力時,誤差主要來源于大氣密度、阻力系數(shù)和迎風(fēng)面積的計(jì)算誤差,阻力系數(shù)和迎風(fēng)面積的誤差可以通過在線辨識的方法予以部分修正[2]。大氣密度主要通過大氣模型計(jì)算得到,但由于影響熱層大氣密度的因素較多且變化機(jī)理復(fù)雜,現(xiàn)有模型都不能準(zhǔn)確描述大氣密度的實(shí)際變化情況,從而造成軌道預(yù)報(bào)誤差。

劉亞英等探索了大氣阻力攝動的精密計(jì)算方法,提高了低軌道衛(wèi)星軌道預(yù)報(bào)精度[3]。Kong等利用三種大氣密度模型,對基于動力學(xué)模型的軌道預(yù)報(bào)方法進(jìn)行了研究[4]。Chen等研究了僅使用兩行軌道根數(shù)提高低軌航天器軌道預(yù)報(bào)精度的方法,利用實(shí)測軌道數(shù)據(jù)對經(jīng)驗(yàn)?zāi)Ｐ瓦M(jìn)行校準(zhǔn)[5]。李夢奇對低軌航天器的大氣參數(shù)修正方法進(jìn)行了研究,降低了大氣阻力攝動誤差并提高了軌道預(yù)報(bào)精度[6]。Ray等利用傅立葉阻力系數(shù)模型對大氣阻力系數(shù)進(jìn)行了精確建模,讓阻力系數(shù)成為與航天器軌道和姿態(tài)相關(guān)參量的函數(shù)[7]。劉舒蒔等基于線性回歸分析建立了大氣阻力系數(shù)補(bǔ)償算法,提高了近地航天器軌道預(yù)報(bào)精度[8]。

強(qiáng)化學(xué)習(xí)是一種機(jī)器學(xué)習(xí)方法,將學(xué)習(xí)看成一個試探-評價的過程[9]。在強(qiáng)化學(xué)習(xí)中,智能體從環(huán)境獲取信息,依據(jù)動作策略采取動作作用于環(huán)境;環(huán)境接受智能體的動作并發(fā)生變化,同時反饋給智能體獎勵或懲罰,即環(huán)境對智能體的動作進(jìn)行評價;智能體根據(jù)環(huán)境反饋信號和當(dāng)前觀測的環(huán)境狀態(tài)選擇下一時刻的動作,依據(jù)是能夠最大化智能體得到的獎勵或者累積獎勵[10]。

諸多學(xué)者對強(qiáng)化學(xué)習(xí)技術(shù)在航天動力學(xué)領(lǐng)域中的應(yīng)用開展了研究。Hovell等提出了一種基于深度強(qiáng)化學(xué)習(xí)的航天器接近操作制導(dǎo)策略,設(shè)計(jì)了一個航天器姿態(tài)跟蹤和對接場景,并驗(yàn)證了該策略的可行性[11]。吳其昌基于MADDPG強(qiáng)化學(xué)習(xí)算法,對航天器追逃博弈問題進(jìn)行了研究,并能滿足航天器機(jī)動策略求解的實(shí)時性與快速性的要求[12]。王月嬌等將DQN強(qiáng)化學(xué)習(xí)算法用于衛(wèi)星自主姿態(tài)控制,有效解決了傳統(tǒng)PD控制器依賴被控對象質(zhì)量參數(shù)的問題[13]。Nicholas等采用強(qiáng)化學(xué)習(xí)算法來生成一種新型閉環(huán)控制器,該控制器能適用于復(fù)雜空間動態(tài)區(qū)域的星載低推力制導(dǎo),并能夠在初始偏差較大的情況下直接引導(dǎo)航天器,提高了制導(dǎo)精度[14]。

綜上所述,對于航天器軌道動力學(xué)模型中大氣阻力攝動建模不準(zhǔn)確的問題,當(dāng)前的研究工作主要是通過解析的方法進(jìn)行參數(shù)矯正,進(jìn)而提高航天器軌道預(yù)報(bào)精度。強(qiáng)化學(xué)習(xí)算法可以通過智能體在線學(xué)習(xí)適應(yīng)動力學(xué)環(huán)境的變化。本文提出了一種基于軌道動力學(xué)模型修正的軌道預(yù)報(bào)方案,該方案將PPO2強(qiáng)化學(xué)習(xí)算法與大氣參數(shù)可調(diào)的動力學(xué)模型相結(jié)合,在動力學(xué)積分預(yù)報(bào)軌道的過程中實(shí)時修正大氣阻力引起的建模誤差,從而提高軌道預(yù)報(bào)精度和預(yù)報(bào)效率。該方案對低軌航天器的軌道預(yù)報(bào)具有一定的參考價值。

2 動力學(xué)模型修正的軌道預(yù)報(bào)方案設(shè)計(jì)

空間站的軌道動力學(xué)方程可表示為

(1)

(2)

式中:Cd為阻力系數(shù);ρ為大氣密度;s為迎風(fēng)面積;m為空間站質(zhì)量;va為空間站相對于大氣的速度。定義空間站的面質(zhì)比與大氣阻力系數(shù)的乘積為彈道系數(shù)B,即

(3)

式中:阻力系數(shù)Cd和面質(zhì)比s/m由地面實(shí)驗(yàn)數(shù)據(jù)及工程經(jīng)驗(yàn)給定。假設(shè)大氣與地球一起旋轉(zhuǎn),則由空間站在地心慣性坐標(biāo)系中的位置r和速度v可得

va=v-?E×r

(4)

式中:?E為地球自轉(zhuǎn)角速度。

軌道預(yù)報(bào)時一般用經(jīng)驗(yàn)大氣模型計(jì)算大氣密度,常用的有US1976標(biāo)準(zhǔn)大氣、NRLMSISE-00模型等。US1976標(biāo)準(zhǔn)大氣模型表示了中等太陽活動期間,由地面到1000km的理性化、靜態(tài)的中緯度平均大氣結(jié)構(gòu),大氣密度的計(jì)算公式為

(5)

式中:Av為阿伏伽德羅常數(shù);mi為N2,O,O2,Ar,He,H六種氣體的標(biāo)稱分子量;Ni為氣體分子的數(shù)密度。Ni只與高度有關(guān),可以通過對標(biāo)準(zhǔn)大氣氣體成分?jǐn)?shù)密度隨高度變化表插值得到。因此基于US1976模型計(jì)算得到的大氣密度只和高度有關(guān),故計(jì)算時效率較高。相比于US1976模型,NRLMSISE-00大氣模型計(jì)算精度較高,其詳細(xì)計(jì)算公式可參考文獻(xiàn)[15]。NRLMSISE-00大氣模型的大氣密度計(jì)算如下:

(6)

在計(jì)算大氣密度時,NRLMSISE-00大氣模型在US1976模型的基礎(chǔ)上添加了6種氣體分子數(shù)密度的修正項(xiàng)ΔNi,并引入了氮N和電離層正氧離子O+兩項(xiàng)?；谥袊臻g站的軌道根數(shù),分別采用兩種大氣模型計(jì)算得到1天內(nèi)空間站所在高度的大氣密度,結(jié)果如圖1所示。

圖1 大氣密度與大氣密度之差隨時間變化

由圖1可知,US1976模型計(jì)算得到的大氣密度隨時間變化幅度較小,NRLMSISE-00模型計(jì)算得到的大氣密度隨時間呈周期性波動,波動周期大致和空間站軌道周期相同。

在實(shí)際工程應(yīng)用中,根據(jù)星載軌道大氣環(huán)境探測器的探測結(jié)果發(fā)現(xiàn),在太陽和地磁活動平靜期,NRLMSISE-00大氣模型計(jì)算得到的大氣密度日均誤差值也能達(dá)到10%,其變化規(guī)律難以用解析表達(dá)式描述。且由于空間站的真實(shí)軌道數(shù)據(jù)難以實(shí)時獲取,故考慮在各時刻基于經(jīng)驗(yàn)大氣模型計(jì)算得到的大氣密度的基礎(chǔ)上引入一個連續(xù)隨機(jī)噪聲,作為真實(shí)大氣密度值,然后進(jìn)行軌道數(shù)值積分,將得到的數(shù)據(jù)作為真實(shí)軌道數(shù)據(jù)。為驗(yàn)證強(qiáng)化學(xué)習(xí)方法的修正效果,本文工作中引入了一個隨機(jī)信號Δk(t)來模擬真實(shí)大氣密度,即

(7)

式中:ρ0為基于NRLMSISE-00大氣模型計(jì)算得到的大氣密度。|Δk(t)|<0.15。真實(shí)大氣密度具有時空連續(xù)性,故Δk(t)是一個隨時間連續(xù)變化的量,考慮采用隨時間呈正態(tài)分布的變量X(t)來描述Δk(t)的變化。X(t)經(jīng)過標(biāo)準(zhǔn)化處理(仿真時段長度為TN,單位為天;規(guī)范化系數(shù)w=0.15),限制Δk(t)的取值范圍。

相比于NRLMSISE-00模型,US1976模型的計(jì)算量更小,因此考慮在軌道預(yù)報(bào)中,通過對US1976模型進(jìn)行修正來計(jì)算大氣密度。根據(jù)式(7),定義

(8)

式中:ΔZ(t)表示US1976模型與真實(shí)大氣模型密度之差。若能通過強(qiáng)化學(xué)習(xí)得到ΔZ(t)隨時間的變化規(guī)律,并補(bǔ)償?shù)筋A(yù)報(bào)模型中,則能有效提高軌道預(yù)報(bào)精度和效率。

基于上述分析,將空間站軌道預(yù)報(bào)的運(yùn)動模型表示為

(9)

式中:ad,agrav,asm,asrp分別為大氣阻力攝動、地球非球形引力攝動、日月引力攝動、太陽光壓攝動產(chǎn)生的攝動加速度,

(10)

記Δβ=BΔZ。方程(9)即為大氣參數(shù)可調(diào)的動力學(xué)模型。Δβ綜合考慮了大氣模型計(jì)算不準(zhǔn)確和軌道動力學(xué)模型建模時彈道系數(shù)B估計(jì)不準(zhǔn)確的影響。由之前的分析可知ΔZ(t)隨時間呈周期性變化,實(shí)際工程應(yīng)用中彈道系數(shù)B的估計(jì)誤差較小,故Δβ隨時間呈周期性變化。若能夠利用強(qiáng)化學(xué)習(xí)網(wǎng)絡(luò)學(xué)習(xí)得到大氣參數(shù)修正量Δβ隨時間的變化規(guī)律,則能夠利用強(qiáng)化學(xué)習(xí)網(wǎng)絡(luò)提高空間站軌道動力學(xué)模型建模精度,從而提高軌道預(yù)報(bào)精度。

由此,設(shè)計(jì)一種基于軌道動力學(xué)模型修正的軌道預(yù)報(bào)方案,其原理如圖2、圖3所示。其中箭頭表示利用獎勵函數(shù)Re(t)對強(qiáng)化學(xué)習(xí)網(wǎng)絡(luò)進(jìn)行訓(xùn)練,將大氣參數(shù)修正量Δβ加到大氣參數(shù)可調(diào)整的動力學(xué)模型上。

圖2 訓(xùn)練段軌道動力學(xué)模型修正方案框架

圖3 預(yù)報(bào)段軌道動力學(xué)模型修正方案框架

圖2表示的是強(qiáng)化學(xué)習(xí)網(wǎng)絡(luò)的訓(xùn)練段。已知J2000.0慣性坐標(biāo)系下t-1時刻空間站的狀態(tài)量X(t-1),根據(jù)式(9)得到下一時刻的狀態(tài)量X(t),并與已知的標(biāo)準(zhǔn)軌道數(shù)據(jù)Xlab(t)作差得到狀態(tài)偏差量ΔX(t),

ΔX(t)=X(t)-Xlab(t)

(11)

圖2中RL-NN為強(qiáng)化學(xué)習(xí)網(wǎng)絡(luò),其輸入為狀態(tài)量X(t)和時間t,輸出為大氣模型參數(shù)調(diào)整量Δβ,其映射關(guān)系為

(12)

為補(bǔ)償大氣模型不準(zhǔn)造成軌道預(yù)報(bào)偏差,應(yīng)使?fàn)顟B(tài)偏差量ΔX接近于0,因此將RL-NN的獎勵函數(shù)R定義為ΔX的函數(shù):

R(t)=f(X(t),Xlab)

(13)

按照圖2完成網(wǎng)絡(luò)訓(xùn)練后,利用強(qiáng)化學(xué)習(xí)網(wǎng)絡(luò)進(jìn)行軌道預(yù)報(bào),其原理如圖3所示。

3 強(qiáng)化學(xué)習(xí)算法設(shè)計(jì)與樣本生成

3.1 強(qiáng)化學(xué)習(xí)過程建模與訓(xùn)練流程設(shè)計(jì)

強(qiáng)化學(xué)習(xí)過程可抽象為馬爾可夫決策過程。其中,狀態(tài)空間S中的狀態(tài)量為[X(t),t]T,動作空間A中動作量為[Δβ],狀態(tài)轉(zhuǎn)移函數(shù)P由大氣參數(shù)可調(diào)的動力學(xué)模型決定;獎勵函數(shù)Re為f(X(t),Xlab),表征t時刻軌道數(shù)值積分結(jié)果X(t)與標(biāo)準(zhǔn)狀態(tài)量Xlab的接近程度,獎勵函數(shù)值越大,表明兩者越接近,即當(dāng)前智能體所采取的策略表現(xiàn)越好。圖3中RL-NN的輸出量Δβ為連續(xù)型變量,故考慮選擇適用于連續(xù)動作空間的高斯策略作為探索策略進(jìn)行策略優(yōu)化。

近端策略優(yōu)化算法,簡稱PPO算法,該算法的特點(diǎn)是實(shí)現(xiàn)復(fù)雜度較低,超參數(shù)調(diào)試較為容易。在PPO算法中智能體每一步訓(xùn)練時都會選擇新策略,通過自適應(yīng)的方式選擇超參數(shù),確保了智能體與環(huán)境交互過程中每一步獲得的獎勵值單調(diào)不減,從而持續(xù)獲得更優(yōu)的策略。

PPO算法包含策略網(wǎng)絡(luò)與評價網(wǎng)絡(luò),采用重要性采樣方法實(shí)現(xiàn)參數(shù)更新。策略網(wǎng)絡(luò)輸出一個均值μ和標(biāo)準(zhǔn)差σ,得到一個由該均值和方差確定的正態(tài)分布。智能體輸出的動作,即大氣模型參數(shù)調(diào)整量Δβ通過對該正態(tài)分布采樣得到,有

(14)

式中:n為動作矢量的維數(shù);Sample函數(shù)為采樣函數(shù)。在PPO網(wǎng)絡(luò)權(quán)重參數(shù)更新過程中,引入了重要性采樣方法,采樣公式如下:

(15)

(16)

從而將f(x)對于策略分布函數(shù)p的期望值轉(zhuǎn)化為相對于策略分布函數(shù)q的期望值。PPO算法中的策略網(wǎng)絡(luò)可等價為策略分布函數(shù),PPO算法中的策略網(wǎng)絡(luò)包含行為策略網(wǎng)絡(luò)和目標(biāo)策略網(wǎng)絡(luò)兩部分,行為策略網(wǎng)絡(luò)輸出的數(shù)據(jù)用于更新目標(biāo)策略網(wǎng)絡(luò)。

為定量評估PPO算法中策略網(wǎng)絡(luò)的性能表現(xiàn),要定義目標(biāo)函數(shù)。先定義優(yōu)勢函數(shù)Aθ′(st,at),其大小與各時刻智能體從環(huán)境中獲得的獎勵值rt相關(guān),反映了當(dāng)前策略得到的大氣參數(shù)補(bǔ)償量Δβ對US1976模型的有效補(bǔ)償程度,即動力學(xué)模型軌道數(shù)值積分得到的狀態(tài)量與標(biāo)準(zhǔn)狀態(tài)量的偏差量越小,則補(bǔ)償效果越好,rt與Aθ′(st,at)的取值越大。PPO2算法是PPO算法的一種,運(yùn)用了截?cái)嗪瘮?shù)來控制策略網(wǎng)絡(luò)參數(shù)的梯度更新,PPO2算法的截?cái)嗪瘮?shù)為

(17)

式中:clip函數(shù)為TensorFlow中numpy函數(shù)庫所包含的函數(shù),用于限制ht(θ)的取值范圍,min為取最小值的函數(shù),θ為策略網(wǎng)絡(luò)參數(shù),式(17)中ht(θ)的定義如下:

(18)

式(17)引入了超參數(shù)ε對重要性采樣獲得的優(yōu)勢函數(shù)進(jìn)行截?cái)嗷幚?在策略網(wǎng)絡(luò)參數(shù)進(jìn)行梯度更新時控制最大更新步長,防止參數(shù)更新速度過快導(dǎo)致訓(xùn)練發(fā)散。PPO2算法進(jìn)行策略網(wǎng)絡(luò)參數(shù)梯度更新的目標(biāo)函數(shù)為

clip(ht(θ),1-ε,1+ε)Aθk(st,at)]

(19)

PPO2算法的網(wǎng)絡(luò)結(jié)構(gòu)與訓(xùn)練流程如圖4所示。PPO2算法仿真訓(xùn)練流程的具體步驟可描述如下。

圖4 PPO2網(wǎng)絡(luò)結(jié)構(gòu)與訓(xùn)練流程

At=R-v

(21)

Step 1:將智能體狀態(tài)量s輸入到策略網(wǎng)絡(luò),得到大氣參數(shù)修正量Δβ的均值μ和標(biāo)準(zhǔn)差σ,由該均值與標(biāo)準(zhǔn)差構(gòu)建一個正態(tài)分布,采樣得到動作量Δβ,將其輸入到動力學(xué)環(huán)境中得到獎勵值rt和下一仿真時刻的狀態(tài)s′,然后存儲該時刻[(s,a,rt),…],將s′輸入到策略網(wǎng)絡(luò),循環(huán)該過程,直到存儲了一定量的[(s,a,rt),…]。

Step 2:將Step1循環(huán)完最后一步得到的s′輸入到評價網(wǎng)絡(luò)中,得到狀態(tài)的v′值,計(jì)算折扣獎勵,即獎勵函數(shù),有

R(t)=rt+γrt+1+…+γTE-t+1rT-1+γTE-tv′

(20)

式中:γ為折扣因子,表示智能體在學(xué)習(xí)過程中對未來潛在獎勵的關(guān)注程度。由此得到各時刻折扣獎勵列表R=[R[0],R[1],…,R[t],…,R[TE]],TE為一個訓(xùn)練仿真回合的最后時刻。

Step 3:將存儲的所有狀態(tài)s組合輸入到評價網(wǎng)絡(luò)中,得到所有狀態(tài)下的v值,計(jì)算優(yōu)勢函數(shù),

Step 4:計(jì)算評價網(wǎng)絡(luò)優(yōu)化的損失函數(shù)c_loss,取均方誤差,

c_loss=MSE(At)

(22)

將c_loss值反向傳播,用于更新評價網(wǎng)絡(luò)。

Step 5:將存儲的所有狀態(tài)s組合輸入策略網(wǎng)絡(luò)1和策略網(wǎng)絡(luò)2,輸出大氣參數(shù)修正量Δβ的均值和標(biāo)準(zhǔn)差,分別構(gòu)造得到正態(tài)分布1和2。將樣本池中所有動作量Δβj組合為動作序列[Δβ1,Δβ2,…]傳入到正態(tài)分布1和2,得到相應(yīng)狀態(tài)下每個Δβj對應(yīng)的概率值prob1和prob2,用prob2除以prob1得到重要性權(quán)重IW。

Step 6:根據(jù)式(19)計(jì)算策略網(wǎng)絡(luò)的損失函數(shù)a_loss,如圖4所示,進(jìn)行反向傳播,更新策略網(wǎng)絡(luò)2。

Step 7:循環(huán)Step 5～6,一定步數(shù)后循環(huán)結(jié)束,用策略網(wǎng)絡(luò)2權(quán)重更新策略網(wǎng)絡(luò)1。

Step 8:循環(huán)Step 1～7,直到滿足迭代結(jié)束條件,即達(dá)到訓(xùn)練精度要求或最大訓(xùn)練輪次,訓(xùn)練結(jié)束。

3.2 強(qiáng)化學(xué)習(xí)訓(xùn)練參量與網(wǎng)絡(luò)模型設(shè)計(jì)

在軌道預(yù)報(bào)場景中,強(qiáng)化學(xué)習(xí)智能體訓(xùn)練時觀察到環(huán)境信息o為大氣參數(shù)可調(diào)的軌道動力學(xué)模型在數(shù)值積分后,得到的空間站各時刻運(yùn)動狀態(tài)量S以及時刻量T,記為

o=[S,T]

(23)

智能體狀態(tài)量s通過軌道數(shù)值積分得到,即引入強(qiáng)化學(xué)習(xí)網(wǎng)絡(luò)的軌道動力學(xué)模型進(jìn)行軌道積分得到空間站各時刻的運(yùn)動狀態(tài)量S,記為

s=S

(24)

智能體的動作為PPO2算法中策略網(wǎng)絡(luò)的輸出,即各時刻大氣參數(shù)調(diào)整量Δβ采樣的正態(tài)分布的均值和標(biāo)準(zhǔn)差,

a=[μ,σ]Δβ

(25)

獎勵函數(shù)選擇為預(yù)期結(jié)果與實(shí)際預(yù)測結(jié)果的范數(shù),用于表征兩者的偏離程度。經(jīng)反復(fù)試驗(yàn)后,在仿真中最終選擇了兩種獎勵函數(shù):

(26)

式中:L為已知的標(biāo)準(zhǔn)軌道數(shù)據(jù)在J2000.0系下各時刻的空間站運(yùn)動狀態(tài)量。獎勵函數(shù)取負(fù)值的原因是PPO2算法的優(yōu)化目標(biāo)是最大化與累加獎勵相關(guān)的目標(biāo)函數(shù)項(xiàng)。

PPO2網(wǎng)絡(luò)包含策略網(wǎng)絡(luò)和評價網(wǎng)絡(luò)。策略網(wǎng)絡(luò)包含目標(biāo)網(wǎng)絡(luò)與行為網(wǎng)絡(luò),且假定兩者具有完全相同的結(jié)構(gòu),包含7個輸入量,即各時刻空間站的運(yùn)動狀態(tài)量及對應(yīng)的時間,策略網(wǎng)絡(luò)包含2個輸出量,即各時刻大氣參數(shù)修正量Δβ正態(tài)分布采樣的均值和方差;評價網(wǎng)絡(luò)包含7個輸入量,與策略網(wǎng)絡(luò)的輸入一致,評價網(wǎng)絡(luò)包含1個輸出量,為當(dāng)前的狀態(tài)值函數(shù)v。

按照從少到多的原則,經(jīng)過反復(fù)試驗(yàn),最終選擇隱藏層數(shù)為3、神經(jīng)元個數(shù)為200的策略網(wǎng)絡(luò)模型,以及隱藏層數(shù)為3、神經(jīng)元個數(shù)為220的評價網(wǎng)絡(luò)模型。隱藏層激活函數(shù)選擇SELU函數(shù),輸出層激活函數(shù)選擇tanh函數(shù),其中SELU函數(shù)的公式如下:

(27)

式中:α和λ是訓(xùn)練前預(yù)設(shè)的。網(wǎng)絡(luò)訓(xùn)練優(yōu)化算法采用Adam算法,性能指標(biāo)選擇為回合累加獎勵?；睾侠塾?jì)獎勵rs為一個仿真回合中,各時刻環(huán)境反饋獎勵ri(i=0,1,2,…,TE)的累加,有

rs=r0+r1+r2+…+rTE

(28)

3.3 訓(xùn)練與測試數(shù)據(jù)生成

軌道預(yù)報(bào)前需要用標(biāo)準(zhǔn)軌道數(shù)據(jù)對強(qiáng)化學(xué)習(xí)網(wǎng)絡(luò)進(jìn)行訓(xùn)練。選擇2022年4月25日某時中國空間站的軌道參數(shù)作為初始標(biāo)準(zhǔn)軌道參數(shù),其半長軸為6765.12km,偏心率為0.0012238,軌道傾角為41.4694°,升交點(diǎn)赤經(jīng)為182.3343°,緯度幅角為210.1109°,平近點(diǎn)角為130.8179°。

在NRLMSISE-00大氣模型上加入了一個時長13天的連續(xù)維納噪聲信號|Δk(t)|<0.15模擬真實(shí)大氣密度,然后通過軌道積分得到模擬的真實(shí)軌道數(shù)據(jù)。仿真時長設(shè)置為13天,引力場模型選擇EGM2008(72×72階),光壓模型選擇BERN模型,三體引力僅考慮日月引力,日月星歷選擇DE405,數(shù)值積分方法選擇RKF7(8),最大積分步長設(shè)置為10s,地磁指數(shù)為3,日均F10.7為270,面質(zhì)比選擇為0.02m2/kg,大氣阻力系數(shù)為2.2。以前10天的數(shù)據(jù)作為訓(xùn)練數(shù)據(jù)樣本,后3天的數(shù)據(jù)作為測試數(shù)據(jù)樣本。

由于原始數(shù)據(jù)的量級差別較大,需要在訓(xùn)練前進(jìn)行歸一化處理。假設(shè)某時刻空間站的運(yùn)動狀態(tài)量為[rx,ry,rz,vx,vy,vz],以仿真時段內(nèi)各方向上位置、速度量模的最大值為標(biāo)準(zhǔn),對策略網(wǎng)絡(luò)與評價網(wǎng)絡(luò)的輸入量進(jìn)行歸一化處理。

4 數(shù)值仿真分析

PPO2算法中對策略網(wǎng)絡(luò)和評價網(wǎng)絡(luò)的權(quán)重參數(shù)按照標(biāo)準(zhǔn)正態(tài)分布進(jìn)行初始化后,還需要對PPO2網(wǎng)絡(luò)訓(xùn)練初始化條件進(jìn)行設(shè)置,各項(xiàng)初始化條件設(shè)置如表1所示。

表1 仿真訓(xùn)練參數(shù)設(shè)置

為保證策略網(wǎng)絡(luò)表現(xiàn)性能的單調(diào)提升,記錄每經(jīng)過一個回合(即對應(yīng)10天訓(xùn)練樣本)訓(xùn)練后訓(xùn)練樣本的累計(jì)獎勵(性能指標(biāo),見式(28))。若性能指標(biāo)相比上一輪提升超過1%,則保存此時的神經(jīng)網(wǎng)絡(luò)模型參數(shù);若新一輪訓(xùn)練后性能指標(biāo)變差,則載入當(dāng)前已獲得的最優(yōu)網(wǎng)絡(luò)模型參數(shù)進(jìn)行訓(xùn)練,直到最終模型性能表現(xiàn)不再有明顯提升。當(dāng)完成一個回合訓(xùn)練后,清零累積獎勵,重置動力學(xué)環(huán)境參數(shù),開始下一回合的訓(xùn)練。

經(jīng)過100回合訓(xùn)練后,訓(xùn)練樣本對應(yīng)累積獎勵的變化情況如圖5所示。采用兩種獎勵函數(shù)時,性能指標(biāo)分別提升92.97%、92.84%,累積獎勵分別從-0.20285提升到-0.01427,從-0.213741提升到-0.015303。

圖5 累積獎勵隨訓(xùn)練回合變化(訓(xùn)練段)

利用測試樣本測試強(qiáng)化學(xué)習(xí)網(wǎng)絡(luò)表現(xiàn)時,強(qiáng)化學(xué)習(xí)網(wǎng)絡(luò)不再更新,直接用于軌道預(yù)報(bào),軌道預(yù)報(bào)結(jié)果與標(biāo)準(zhǔn)軌道數(shù)據(jù)之間存在的偏差即強(qiáng)化學(xué)習(xí)方法的軌道預(yù)報(bào)精度,由此驗(yàn)證強(qiáng)化學(xué)習(xí)方法的泛化能力。

基于測試樣本,分別采用兩種獎勵函數(shù)進(jìn)行軌道預(yù)報(bào)時的軌道根數(shù)誤差變化如圖6所示,進(jìn)行軌道預(yù)報(bào)的位置誤差(J2000.0系下)與僅考慮NRLMSISE-00大氣模型進(jìn)行軌道預(yù)報(bào)時的位置誤差以及程序運(yùn)行時間的對比情況如圖7、圖8、表2所示。

表2 基于經(jīng)驗(yàn)大氣模型與混合動力學(xué)模型的最大預(yù)報(bào)位置誤差與預(yù)報(bào)時間

續(xù)圖6

圖7 基于經(jīng)驗(yàn)大氣模型與混合動力學(xué)模型所得位置誤差(獎勵函數(shù)一)

圖8 基于經(jīng)驗(yàn)大氣模型與混合動力學(xué)模型所得位置誤差(獎勵函數(shù)二)

由圖6、圖7、圖8和表2可知,在測試樣本中,相比于僅考慮NRLMSISE-00模型,包含強(qiáng)化學(xué)習(xí)網(wǎng)絡(luò)的混合動力學(xué)模型用于軌道預(yù)報(bào)的誤差更小。采用獎勵函數(shù)一時X、Y、Z方向上的最大位置誤差與最大總位置誤差分別減少了94.21%、92.62%、93.09%、93.19%;在采用獎勵函數(shù)二時X、Y、Z方向上的最大位置誤差與最大總位置誤差分別減少了93.91%、91.96%、92.79%、92.81%,且軌道根數(shù)誤差變化更為平緩,波動幅度較小。將訓(xùn)練后得到的強(qiáng)化學(xué)習(xí)網(wǎng)絡(luò)用于預(yù)報(bào),程序運(yùn)行時間分別減小了83.40%,83.29%,由此驗(yàn)證了強(qiáng)化學(xué)習(xí)方法的有效性。

另構(gòu)造與PPO2網(wǎng)絡(luò)中策略網(wǎng)絡(luò)隱藏層層數(shù)與神經(jīng)元個數(shù)相同的深度神經(jīng)網(wǎng)絡(luò)模型,令其輸入量為僅考慮NRLMSISE-00模型進(jìn)行軌道數(shù)值積分得到的空間站J2000.0系下各時刻的狀態(tài)量,令其輸出量為輸入狀態(tài)量相比于標(biāo)準(zhǔn)狀態(tài)量的偏差量,采用10天的標(biāo)準(zhǔn)軌道數(shù)據(jù)對神經(jīng)網(wǎng)絡(luò)模型進(jìn)行訓(xùn)練,隨后進(jìn)行3天的軌道預(yù)報(bào),預(yù)報(bào)狀態(tài)量為神經(jīng)網(wǎng)絡(luò)輸出狀態(tài)偏差量與神經(jīng)網(wǎng)絡(luò)輸入狀態(tài)量的和,此時軌道預(yù)報(bào)的最大位置誤差與程序運(yùn)行時間如表3所示。

表3 基于強(qiáng)化學(xué)習(xí)網(wǎng)絡(luò)模型的最大預(yù)報(bào)位置誤差與預(yù)報(bào)時間

由表3可知,在進(jìn)行軌道預(yù)報(bào)時,相比于強(qiáng)化學(xué)習(xí)方法,直接利用深度神經(jīng)網(wǎng)絡(luò)對原動力學(xué)模型的軌道數(shù)值積分結(jié)果偏差進(jìn)行補(bǔ)償,其預(yù)報(bào)的程序運(yùn)行時間更長,預(yù)報(bào)精度更低,由此驗(yàn)證了本文研究的強(qiáng)化學(xué)習(xí)方法的優(yōu)越性。

5 結(jié)論

本文針對大氣阻力攝動難以準(zhǔn)確建模的問題,提出了一種融合PPO2強(qiáng)化學(xué)習(xí)網(wǎng)絡(luò)和軌道動力學(xué)模型修正的軌道預(yù)報(bào)方案,隨后完成了網(wǎng)絡(luò)模型構(gòu)建,訓(xùn)練和測試。結(jié)合最終仿真測試結(jié)果進(jìn)行如下總結(jié):

1)該方案具有較高的預(yù)報(bào)精度和預(yù)報(bào)效率,相比于經(jīng)驗(yàn)大氣模型,采用兩種獎勵函數(shù)時,三方向上最大位置誤差與最大總位置誤差分別減少了94.21%、92.62%、93.09%、93.19%及93.91%、91.96%、92.79%、92.81%;進(jìn)行軌道預(yù)報(bào)時,程序運(yùn)行時間分別減小了83.40%、83.29%。

2)利用PPO2算法和強(qiáng)化學(xué)習(xí)網(wǎng)絡(luò)可有效修正經(jīng)驗(yàn)大氣模型得到的大氣模型參數(shù)存在的誤差,與深度學(xué)習(xí)方法不同,強(qiáng)化學(xué)習(xí)方法引入了獎勵函數(shù),在與動力學(xué)環(huán)境互動過程中為強(qiáng)化學(xué)習(xí)網(wǎng)絡(luò)的訓(xùn)練優(yōu)化提供了一個學(xué)習(xí)目標(biāo)。

3)與空間站軌道預(yù)報(bào)類似,該方案對其他低軌航天器的軌道預(yù)報(bào)同樣具有一定參考價值。

未來的研究工作中,可以考慮在本文研究工作的基礎(chǔ)上,利用高精度定軌得到的軌道數(shù)據(jù)用于大氣模型參數(shù)的實(shí)時修正。