基于經驗模式分解和多項式擬合的靜止衛星變軌預測

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

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

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

常用的噪聲濾除方法[6-9]為使用小波變換在不同尺度上分解帶噪信號，再在不含噪尺度上對信號進行分析，但其性能受小波基選取影響。經驗模式分解方法[10]可根據信號本身的時間尺度特性自適應選取“基函數”，克服了小波基選擇難題。

該文在經驗模式分解域，通過殘差分量多項式預測值與真實值誤差的閾值化處理，實時判斷靜止衛星是否變軌。

1 EMD方法和一致性EMD方法

經驗模式分解（Empirical Mode Decomposition，EMD）方法能夠根據信號本身的特性，自適應地將時域信號分解為一組不同尺度的固有模態函數（Intrinsic Mode Function，IMF），而每個尺度上分解所得殘差分量代表原始信號的趨勢或均值。

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

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

2 該文變軌預測算法實施步驟

經驗模式分解和多項式擬合的靜止衛星變軌預測算法的基本思想是利用EMD 方法的多尺度濾波特性[12-14]，對衛星經度信號在多個尺度上進行平滑濾波，并對濾波后的信號進行多項式擬合、預測，最終通過預測誤差判斷衛星是否變軌。

2.1 經度信號預處理

在靜止軌道衛星發射階段，其經度信號會產生大幅度變化。這種劇烈變化會導致多項式擬合誤差增大，影響靜止衛星變軌判斷算法的準確度。圖1（a）顯示了某靜止衛星的歷史經度信號。圖中衛星發射階段經度信號、運行周期信號同步變化，且變化幅度遠大于其進入靜止軌道后正常變軌時變化幅度。如果單一采用靜止衛星經度信號進行變軌預測，會造成將發射階段的經度變化誤檢測為正常變軌情況。

圖1 衛星經度信號預處理效果

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

2.2 預處理后經度信號經驗模式分解

經驗模式分解方法得到的各階IMF 分量cn(i)包含了待分解信號x(i)或rn-1(i)中高頻成分，rn(i)則包含了x(i)或rn-1(i)中低頻成分，因此rn(i)可看作原始信號x(i)經不同平滑濾波器輸出的、在不同尺度上的數據。階數n越小則rn(i)就越接近原始信號x(i)，所含噪聲成分越多，變軌預測算法受噪聲影響越大；階數n越大則rn(i)就越平滑，所含噪聲成分越少，變軌預測算法受噪聲影響越小；同時，階數n過大、rn(i)過分平滑會導致rn(i)丟失細節而導致變軌預測結果錯誤，因此合理選擇后續擬合分析信號rn(i)的階數n，直接影響最終變軌預測準確性。該文選取第1 尺度上的殘差分量r1(i)作為擬合分析信號。同時，由于該算法是預測算法，經驗模式分解、多項式擬合預測等全部操作均基于1 到i-1 時刻的歷史數值，因此每當獲取到新的x(i)值后，需對x(1:i)進行經驗模式分解，得到r1(i)；r1(i)并非由x(i)一次經驗模式分解計算得出。

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

圖2 預處理后經度信號的分解結果

2.3 多項式擬合與預測

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

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

2.4 噪聲閾值估計與變軌判決

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

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

3 實驗結果與分析

該文使用Matlab 軟件評估經驗模式分解和多項式擬合的靜止衛星變軌預測算法。其中，EMD 方法的篩選門限SD 取0.2；圖1（a）所示靜止衛星的經度信號作為待分析信號；圖1（c）所示預處理后經度信號x(i)的兩個變軌事件分別發生在第81-125、277-297 點之間。

圖3 給出了預處理后經度信號x(i)采用經驗模式分解和多項式擬合的靜止衛星變軌預測算法預測變軌的例子。其中，圖3（a）顯示的是靜止衛星在第81 個時間點開始變軌的預測結果，圖3（b）顯示的是靜止衛星在第277 個時間點開始變軌的預測結果。圖中，星號曲線表示第1 尺度上的殘差分量實際值（曲線中9 個點取值均由預處理后經度信號x(1∶81)或x(1∶277)采用EMD 方法計算得出），空心圓曲線表示第1 尺度上的殘差分量預測值（曲線中前8 個點取值由預處理后經度信號x(1∶80)或x(1∶276)采用EMD方法計算得出，第9 個點取值由前8 個點的值擬合預測得出，即re1(81)或re1(277)），實心圓曲線表示第1 尺度上的變軌判決門限（曲線中9 個點取值均由re1(81)或re1(277)減去噪聲幅度門限Tj得出）。由圖可見，在第81、227 這兩個采樣點，第1 尺度上的殘差分量實際值超過了變軌判決門限，該算法對靜止軌道衛星變軌事件實現了準確、及時預測。

圖3 該文算法的變軌預測結果

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

圖4 軌位小幅變化的預測結果比對

4 結束語

該文提出一種經驗模式分解和多項式擬合的靜止衛星變軌預測算法。該算法基于EMD 方法的多尺度濾波特性，將殘差信號的多項式預測誤差與噪聲閾值對比，實現從衛星經度信號中預測衛星變軌事件。以實際靜止衛星為例的實驗表明，該算法能夠有效對靜止衛星變軌進行實時預測，為實現變軌事件的自動預測指明了方向。