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

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

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

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

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

1 EMD方法和一致性EMD方法

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

設(shè)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 方法具有尺度濾波特性。

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

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

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

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

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

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

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

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

經(jīng)驗?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ù)測算法受噪聲影響越?。煌瑫r，階數(shù)n過大、rn(i)過分平滑會導致rn(i)丟失細節(jié)而導致變軌預(yù)測結(jié)果錯誤，因此合理選擇后續(xù)擬合分析信號rn(i)的階數(shù)n，直接影響最終變軌預(yù)測準確性。該文選取第1 尺度上的殘差分量r1(i)作為擬合分析信號。同時，由于該算法是預(yù)測算法，經(jīng)驗?zāi)Ｊ椒纸?、多項式擬合預(yù)測等全部操作均基于1 到i-1 時刻的歷史數(shù)值，因此每當獲取到新的x(i)值后，需對x(1:i)進行經(jīng)驗?zāi)Ｊ椒纸?，得到r1(i)；r1(i)并非由x(i)一次經(jīng)驗?zāi)Ｊ椒纸庥嬎愕贸觥?/p>

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

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

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

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

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

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

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

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

3 實驗結(jié)果與分析

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

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

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

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

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

4 結(jié)束語

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