基于經(jīng)驗(yàn)?zāi)Ｊ椒纸夂投囗?xiàng)式擬合的靜止衛(wèi)星變軌預(yù)測

李 燁，石會鵬，張 周，上官澤胤

（1.31007 部隊(duì)，北京 100079；2.北京東方波泰無線電頻譜技術(shù)研究所有限公司，北京 100041；3.軍事科學(xué)院，北京 100071）

空間頻率軌道資源按照先登先占原則申請使用[1]。為維護(hù)我國空間頻率軌道資源安全、拓展空間頻率軌道資源利益，掌握全球靜止衛(wèi)星的軌位數(shù)據(jù)成為關(guān)鍵[2-3]。受姿態(tài)控制、攝動力等因素影響，實(shí)際運(yùn)行中靜止衛(wèi)星的軌位具有小幅波動[4]。從具有波動的經(jīng)度信號等軌道參數(shù)中，及時(shí)、準(zhǔn)確識別并預(yù)測靜止衛(wèi)星變軌，仍面臨缺少自動化、高魯棒性算法的困難[5]。

常用的噪聲濾除方法[6-9]為使用小波變換在不同尺度上分解帶噪信號，再在不含噪尺度上對信號進(jìn)行分析，但其性能受小波基選取影響。經(jīng)驗(yàn)?zāi)Ｊ椒纸夥椒╗10]可根據(jù)信號本身的時(shí)間尺度特性自適應(yīng)選取“基函數(shù)”，克服了小波基選擇難題。

該文在經(jīng)驗(yàn)?zāi)Ｊ椒纸庥?，通過殘差分量多項(xiàng)式預(yù)測值與真實(shí)值誤差的閾值化處理，實(shí)時(shí)判斷靜止衛(wèi)星是否變軌。

1 EMD方法和一致性EMD方法

經(jīng)驗(yàn)?zāi)Ｊ椒纸猓‥mpirical Mode Decomposition，EMD）方法能夠根據(jù)信號本身的特性，自適應(yīng)地將時(shí)域信號分解為一組不同尺度的固有模態(tài)函數(shù)（Intrinsic Mode Function，IMF），而每個(gè)尺度上分解所得殘差分量代表原始信號的趨勢或均值。

設(shè)x(i)為待分解的原始信號，i=1,…,I；cn(i)為第n次分解得到的IMF 分量；rn(i)為第n次分解得到的殘差分量，n=1，…,N。EMD 算法將原始信號x(i)分解為N個(gè)IMF 分量c1(i)、c2(i)、…、cN(i)和一個(gè)殘差分量rN(i)的和。其中，N個(gè)IMF 分量分別包含了原信號x(i)中不同時(shí)間特征尺度大小的成分，且各IMF分量的時(shí)間特征尺度隨n增加而依次變大，這說明EMD 方法具有尺度濾波特性。

在邊緣提取算法中，一般選擇暗區(qū)域或噪聲信號在某個(gè)尺度上信號的特征量作為該尺度上檢測閾值。但是，如果分別對待分析信號和噪聲進(jìn)行經(jīng)驗(yàn)?zāi)Ｊ椒纸?，由于EMD 方法是自適應(yīng)基函數(shù)分解方法，且兩信號的尺度特性不同，那么就不能保證兩個(gè)信號的第n階IMF 分量在同一尺度上，也就不能使用噪聲信號第n階IMF 或殘差分量的特征量作為該尺度上邊緣檢測的自動化判決閾值。文獻(xiàn)[11]利用EMD 方法中計(jì)算上、下包絡(luò)線的內(nèi)插操作僅與擬合點(diǎn)鄰近有限個(gè)極值點(diǎn)（例如三次樣條曲線內(nèi)插需要臨近5 個(gè)極大值點(diǎn)或5 個(gè)極小值點(diǎn)）有關(guān)的特點(diǎn)，通過保證參考噪聲各階IMF 分量的篩分次數(shù)與原始信號x(i)各階IMF 分量的篩分次數(shù)相同的方法，來估計(jì)各階rn(i)中噪聲信號。這種對參考噪聲的分解方法既保證了信號x(i)與參考噪聲分解得到的IMF 分量數(shù)相同，又保證了兩個(gè)信號的第n階IMF 分量在同一個(gè)尺度上，即保證了兩者的一致性，因此稱其為一致性經(jīng)驗(yàn)?zāi)Ｊ椒纸猓╟onsistent EMD）方法。

2 該文變軌預(yù)測算法實(shí)施步驟

經(jīng)驗(yàn)?zāi)Ｊ椒纸夂投囗?xiàng)式擬合的靜止衛(wèi)星變軌預(yù)測算法的基本思想是利用EMD 方法的多尺度濾波特性[12-14]，對衛(wèi)星經(jīng)度信號在多個(gè)尺度上進(jìn)行平滑濾波，并對濾波后的信號進(jìn)行多項(xiàng)式擬合、預(yù)測，最終通過預(yù)測誤差判斷衛(wèi)星是否變軌。

2.1 經(jīng)度信號預(yù)處理

在靜止軌道衛(wèi)星發(fā)射階段，其經(jīng)度信號會產(chǎn)生大幅度變化。這種劇烈變化會導(dǎo)致多項(xiàng)式擬合誤差增大，影響靜止衛(wèi)星變軌判斷算法的準(zhǔn)確度。圖1（a）顯示了某靜止衛(wèi)星的歷史經(jīng)度信號。圖中衛(wèi)星發(fā)射階段經(jīng)度信號、運(yùn)行周期信號同步變化，且變化幅度遠(yuǎn)大于其進(jìn)入靜止軌道后正常變軌時(shí)變化幅度。如果單一采用靜止衛(wèi)星經(jīng)度信號進(jìn)行變軌預(yù)測，會造成將發(fā)射階段的經(jīng)度變化誤檢測為正常變軌情況。

圖1 衛(wèi)星經(jīng)度信號預(yù)處理效果

為了消除靜止軌道衛(wèi)星發(fā)射階段對變軌預(yù)測算法影響，該文采用衛(wèi)星運(yùn)行周期信號作為參考信號，以運(yùn)行周期大于23 小時(shí)56 分為判斷依據(jù)，在變軌預(yù)測前對原始經(jīng)度信號進(jìn)行預(yù)處理，將靜止衛(wèi)星發(fā)射階段的經(jīng)度信號從原始經(jīng)度信號中刪除。圖1（c）顯示了某靜止衛(wèi)星經(jīng)度信號預(yù)處理的結(jié)果。

2.2 預(yù)處理后經(jīng)度信號經(jīng)驗(yàn)?zāi)Ｊ椒纸?/h3>

經(jīng)驗(yàn)?zāi)Ｊ椒纸夥椒ǖ玫降母麟AIMF 分量cn(i)包含了待分解信號x(i)或rn-1(i)中高頻成分，rn(i)則包含了x(i)或rn-1(i)中低頻成分，因此rn(i)可看作原始信號x(i)經(jīng)不同平滑濾波器輸出的、在不同尺度上的數(shù)據(jù)。階數(shù)n越小則rn(i)就越接近原始信號x(i)，所含噪聲成分越多，變軌預(yù)測算法受噪聲影響越大；階數(shù)n越大則rn(i)就越平滑，所含噪聲成分越少，變軌預(yù)測算法受噪聲影響越??；同時(shí)，階數(shù)n過大、rn(i)過分平滑會導(dǎo)致rn(i)丟失細(xì)節(jié)而導(dǎo)致變軌預(yù)測結(jié)果錯誤，因此合理選擇后續(xù)擬合分析信號rn(i)的階數(shù)n，直接影響最終變軌預(yù)測準(zhǔn)確性。該文選取第1 尺度上的殘差分量r1(i)作為擬合分析信號。同時(shí)，由于該算法是預(yù)測算法，經(jīng)驗(yàn)?zāi)Ｊ椒纸狻⒍囗?xiàng)式擬合預(yù)測等全部操作均基于1 到i-1 時(shí)刻的歷史數(shù)值，因此每當(dāng)獲取到新的x(i)值后，需對x(1:i)進(jìn)行經(jīng)驗(yàn)?zāi)Ｊ椒纸?，得到r1(i)；r1(i)并非由x(i)一次經(jīng)驗(yàn)?zāi)Ｊ椒纸庥?jì)算得出。

圖2 顯示了對某靜止衛(wèi)星的預(yù)處理后經(jīng)度信號采用EMD 方法進(jìn)行多尺度分解的例子。其中，原始信號x(i)采用EMD 方法分解得到5 個(gè)IMF 分量以及r1(i)至r5(i)共5 個(gè)殘差分量。

圖2 預(yù)處理后經(jīng)度信號的分解結(jié)果

2.3 多項(xiàng)式擬合與預(yù)測

為了實(shí)時(shí)檢測靜止衛(wèi)星變軌，該文采用對比某時(shí)刻某尺度上殘差分量預(yù)測值與真實(shí)值方法實(shí)現(xiàn)。其中，預(yù)測值由某時(shí)刻i+1 之前的歷史殘差分量數(shù)據(jù)，采用滑動均值法求得。參考文獻(xiàn)[15]引述，該文選用0 階多項(xiàng)式對第一尺度上殘差分量r1(i) 中第i-7、i-6、i-5、i-4、i-3、i-2、i-1、i共8 個(gè)時(shí)刻的數(shù)值進(jìn)行擬合，并采用該0 階多項(xiàng)式預(yù)測第i+1時(shí)刻的殘差分量數(shù)值re1(i+1)，即：

與文獻(xiàn)[15]引述方法不同，為了實(shí)現(xiàn)實(shí)時(shí)預(yù)測靜止衛(wèi)星變軌，該文的多項(xiàng)式擬合只使用了后置多項(xiàng)式擬合，避免了前置多項(xiàng)式擬合對某時(shí)刻i+1 之后未來數(shù)據(jù)的依賴。

2.4 噪聲閾值估計(jì)與變軌判決

噪聲幅度門限Tj的選取決定閾值化處理的效果，直接影響變軌預(yù)測算法自動化判決準(zhǔn)確性。如果噪聲幅度門限Tj取值較小，那么無法去除小幅機(jī)動造成的噪聲影響，提高了變軌預(yù)測算法虛警率；如果噪聲幅度門限Tj取值較大，那么將濾除掉部分變軌事件造成的小幅預(yù)測誤差，提高了變軌預(yù)測算法漏警率。參考文獻(xiàn)[16-17]所述，小波域?yàn)V波選取圖像信號中“暗”區(qū)域的信號作為參考噪聲,并通過對參考噪聲進(jìn)行小波變換求得各尺度上噪聲閾值。但由于EMD 是一種自適應(yīng)基函數(shù)濾波方法,單獨(dú)對參考噪聲直接采用EMD 方法分解既不能保證原始信號x(i)與參考噪聲分解得到的分量數(shù)量相同,又不能保證原始信號x(i)和參考噪聲的第n階分量在同一尺度上。因此該文EMD 域變軌判決的噪聲閾值選取方法不能簡單借鑒小波域?yàn)V波。該文首先選取靜止衛(wèi)星在某軌位駐留期間的、平穩(wěn)的經(jīng)度信號作為參考噪聲源信號，并對其去均值處理得到參考噪聲信號；然后采用一致性經(jīng)驗(yàn)?zāi)Ｊ椒纸夥椒▽⒖荚肼曔M(jìn)行分解，得到各尺度上參考噪聲信號,從而計(jì)算出各尺度上參考噪聲信號的功率；最后參考噪聲功率的平方根即為該尺度上噪聲幅度門限Tj。

判決靜止衛(wèi)星當(dāng)前是否正在進(jìn)行變軌的判據(jù)：如果第i+1 時(shí)刻的殘差分量預(yù)測值re1(i+1)與第i+1時(shí)刻的殘差分量真實(shí)值r1(i+1)的差值絕對值不大于第1 尺度上噪聲幅度門限Tj，則判斷在第i+1 時(shí)刻衛(wèi)星未實(shí)施變軌；如果第i+1 時(shí)刻的殘差分量預(yù)測值re1(i+1)與第i+1 時(shí)刻的殘差分量真實(shí)值r1(i+1)的差值絕對值大于第1 尺度上噪聲幅度門限Tj，則判斷在第i+1 時(shí)刻衛(wèi)星可能正在實(shí)施變軌。

3 實(shí)驗(yàn)結(jié)果與分析

該文使用Matlab 軟件評估經(jīng)驗(yàn)?zāi)Ｊ椒纸夂投囗?xiàng)式擬合的靜止衛(wèi)星變軌預(yù)測算法。其中，EMD 方法的篩選門限SD 取0.2；圖1（a）所示靜止衛(wèi)星的經(jīng)度信號作為待分析信號；圖1（c）所示預(yù)處理后經(jīng)度信號x(i)的兩個(gè)變軌事件分別發(fā)生在第81-125、277-297 點(diǎn)之間。

圖3 給出了預(yù)處理后經(jīng)度信號x(i)采用經(jīng)驗(yàn)?zāi)Ｊ椒纸夂投囗?xiàng)式擬合的靜止衛(wèi)星變軌預(yù)測算法預(yù)測變軌的例子。其中，圖3（a）顯示的是靜止衛(wèi)星在第81 個(gè)時(shí)間點(diǎn)開始變軌的預(yù)測結(jié)果，圖3（b）顯示的是靜止衛(wèi)星在第277 個(gè)時(shí)間點(diǎn)開始變軌的預(yù)測結(jié)果。圖中，星號曲線表示第1 尺度上的殘差分量實(shí)際值（曲線中9 個(gè)點(diǎn)取值均由預(yù)處理后經(jīng)度信號x(1∶81)或x(1∶277)采用EMD 方法計(jì)算得出），空心圓曲線表示第1 尺度上的殘差分量預(yù)測值（曲線中前8 個(gè)點(diǎn)取值由預(yù)處理后經(jīng)度信號x(1∶80)或x(1∶276)采用EMD方法計(jì)算得出，第9 個(gè)點(diǎn)取值由前8 個(gè)點(diǎn)的值擬合預(yù)測得出，即re1(81)或re1(277)），實(shí)心圓曲線表示第1 尺度上的變軌判決門限（曲線中9 個(gè)點(diǎn)取值均由re1(81)或re1(277)減去噪聲幅度門限Tj得出）。由圖可見，在第81、227 這兩個(gè)采樣點(diǎn)，第1 尺度上的殘差分量實(shí)際值超過了變軌判決門限，該算法對靜止軌道衛(wèi)星變軌事件實(shí)現(xiàn)了準(zhǔn)確、及時(shí)預(yù)測。

圖3 該文算法的變軌預(yù)測結(jié)果

為了對比該文變軌預(yù)測算法與文獻(xiàn)[15]引述直接對預(yù)處理后經(jīng)度信號x(i)進(jìn)行多項(xiàng)式擬合預(yù)測的算法，實(shí)驗(yàn)進(jìn)一步采用相同的參考噪聲和上述兩種算法，對預(yù)處理后經(jīng)度信號x(i)第273 點(diǎn)的軌位小幅變化進(jìn)行效果分析，結(jié)果如圖4 所示（圖中各曲線代表的內(nèi)容同圖3）。其中，圖4（a）所示為采用該文變軌預(yù)測算法的預(yù)測結(jié)果；圖4（b）所示為直接采用多項(xiàng)式擬合預(yù)測算法的預(yù)測結(jié)果。由圖可見，該文所述變軌預(yù)測算法未在第273 點(diǎn)預(yù)測變軌；直接多項(xiàng)式擬合預(yù)測算法則出現(xiàn)變軌誤判，產(chǎn)生虛警。該文所述變軌預(yù)測算法未出現(xiàn)虛警的原因應(yīng)為EMD 算法消除了小幅軌位變化對多項(xiàng)式擬合預(yù)測的影響。

圖4 軌位小幅變化的預(yù)測結(jié)果比對

4 結(jié)束語

該文提出一種經(jīng)驗(yàn)?zāi)Ｊ椒纸夂投囗?xiàng)式擬合的靜止衛(wèi)星變軌預(yù)測算法。該算法基于EMD 方法的多尺度濾波特性，將殘差信號的多項(xiàng)式預(yù)測誤差與噪聲閾值對比，實(shí)現(xiàn)從衛(wèi)星經(jīng)度信號中預(yù)測衛(wèi)星變軌事件。以實(shí)際靜止衛(wèi)星為例的實(shí)驗(yàn)表明，該算法能夠有效對靜止衛(wèi)星變軌進(jìn)行實(shí)時(shí)預(yù)測，為實(shí)現(xiàn)變軌事件的自動預(yù)測指明了方向。