“北斗”導航衛星定軌殘差特征提取與分布檢驗

陳略 唐歌實 崔紅正 陳明 劉薈萃 王美

(1 北京航天飛行控制中心，北京 100094) (2 航天飛行動力學技術國家級重點實驗室，北京 100094)

“北斗”導航衛星定軌殘差特征提取與分布檢驗

陳略1,2唐歌實1,2崔紅正1,2陳明1,2劉薈萃1,2王美1,2

(1 北京航天飛行控制中心，北京 100094) (2 航天飛行動力學技術國家級重點實驗室，北京 100094)

針對“北斗”導航衛星定軌殘差中非建模系統誤差問題，提出基于總體平均經驗模式分解(EEMD)與Hilbert譜分析相結合的定軌殘差特征提取方法，并對定軌殘差正態分布特性進行卡方檢驗。首先，分析總體經驗模式分解原理，提出濾波輔助的改進EEMD方法與Hilbert譜特征提取結合，建立定軌殘差特征提取模型，并闡述卡方檢驗正態分布原理；然后對北斗導航衛星偽距與載波相位數據定軌殘差特征進行分析；最后將提出的方法應用于多測站、多北斗GEO衛星的定軌殘差分析中。結果表明，EEMD與Hilbert譜方法有效提取出了定軌殘差中1天的軌道周期項，且EEMD處理后定軌殘差的卡方統計量為5.5，其值小于卡方檢驗臨界，可視為正態分布。該方法可為北斗定軌殘差中非建模系統誤差分離、提高定軌內符合精度提供有效技術手段。

定軌殘差；特征提??；分布檢驗；偽距；載波相位；總體平均經驗模式分解；北斗導航衛星系統

1 引言

“北斗”衛星的觀測量包括偽距觀測量與載波相位觀測量，利用偽距觀測量、載波相位觀測量可對“北斗”導航衛星進行高精度定軌。定軌殘差是定軌內符合精度表現形式之一，通過檢測定軌殘差可評估定軌精度[1]。提高衛星星歷、鐘差及相關模型參數的精度，進而提高衛星精密定軌精度，是提高北斗衛星導航系統精密定位服務能力的關鍵[2]。影響北斗定軌解算的非建模系統誤差產生的原因很多，包括多路徑影響、衛星軌道、電離層延遲、對流層延遲及衛星鐘差等。通常電離層、對流層、衛星鐘差等誤差可通過建模進行誤差修正，但對于某些難以模型化的誤差則通常簡單地忽略其影響，這些非建模系統誤差成為高精度北斗定軌的主要誤差源之一，且在定軌殘差中普遍存在。如果能對“北斗”衛星的定軌殘差進行特征提取，找出其隱藏信息，對于“北斗”定軌殘差中非建模系統誤差分離與建模具有重要意義。

本文研究將總體平均經驗模式分解(EEMD)方法引入到定軌殘差分析中，提出濾波輔助的改進EEMD方法，用于精確提取“北斗”定軌殘差中的特征量。EEMD方法來源于經驗模式分解(EMD)[3]，已被迅速成功應用于雷達信號處理、遙感數據處理、生物信號處理、心電信號處理、故障診斷和表面肌電信號處理等領域[4-5]。

針對“北斗”導航衛星定軌殘差分析這一研究問題，本文提出利用EEMD與Hilbert譜的北斗定軌殘差非系統建模誤差特征提取方法，用于精確提取“北斗”定軌殘差中的特征信息，并通過卡方檢驗方法來對定軌殘差的分布特性進行定量檢驗，以期后續能對“北斗”定軌殘差中非建模系統誤差進行準確分離建模，旨在提高對“北斗”導航衛星的精密定軌精度。

2 定軌殘差特征提取

“北斗”定軌殘差是指針對“北斗”導航衛星的觀測值減去觀測量的計算值，定軌殘差直接反應了對衛星定軌的內符合精度。定軌殘差表征是一個時間序列，可看成一個時間序列的信號，因此可以利用時間序列分析方法去分析定軌殘差。而卡方檢驗是檢驗時間序列是否呈正態分布的有效方法。下面介紹“北斗”定軌殘差特征提取與殘差分布檢驗基本理論。

2.1 EEMD分解

EEMD分解的原理步驟[4]如下：首先，通過對信號多次加入具有均勻尺度、幅值標準差為常數的隨機高斯白噪聲序列后，使原信號具有足夠的極值點。然后，對每次加入高斯白噪聲后的信號進行EMD分解得到相應基本模式分量(IMF)，利用不相關隨機序列的統計均值為零的原理，將所有對應的IMF做總體平均運算，消除了多次加入的高斯白噪聲對真實IMF的影響。最后將總體平均后的IMF作為EEMD分解的IMF最終結果，這就是EEMD分解方法[7]。EEMD分解方法有效避免了EMD方法中模式混淆現象。

在EEMD分解方法中，加入的白噪聲(σ)與總體平均次數(N)是兩個最重要參數。加入的高斯白噪聲應既不影響信號中的有效高頻成分極值點間隔的分布特性，同時又能改變信號中低頻成分的極值點間隔的分布特性[5]。文獻[4]提出EEMD方法中加入高斯白噪聲準則：

(1)

圖1 濾波輔助的改進EEMD算法流程Fig.1 Algorithm flow chat of improved EEMD with filter assisting

為更有效地確定σ值，本文提出通過濾波方式來確定σ值，即首先對原始分析信號進行高通或帶通濾波，評估指定頻帶內的信號分布特性，確定εh值，得到σ值，進而確定α值。濾波器通帶頻率選取可根據分析信號特征確定，例如通帶選擇可采用[3fs/8,fs/2]，即包含分析信號的高頻段成分，fs為采樣頻率。這樣就得到一種濾波輔助的改進EEMD方法，其算法流程如圖1所示。

2.2 Hilbert包絡譜

圖2 “北斗”定軌殘差特征提取方法模型Fig.2 Feature extraction flow chart of COMPASS orbit residuals

Hilbert變換的定義[8]為

(2)

基于以上原理，本文提出了總體平均經驗模式分解與Hilbert譜的北斗定軌殘差非系統建模誤差特征提取方法，其處理流程如圖2所示。

2.3 卡方檢驗正態分布原理

卡方統計量計算公式為

(3)

式中Ai為i水平的觀測頻數；Ei為i水平的期望頻數；n為總頻數；pi為i水平的期望頻率；k為單元格數。i水平的期望頻數Ei等于總頻數n×i水平的期望概率pi。

由卡方統計量計算公式可知，當觀察頻數與期望頻數完全一致時，χ2值為0；觀察頻數與期望頻數越接近，兩者之間的差異越小，χ2值越?。环粗?，觀察頻數與期望頻數差別越大，兩者之間的差異越大，χ2值越大。在進行卡方正態分布檢驗時，其方法為：將計算χ2值與理論χ2臨界值進行比較，如果計算χ2值小于理論χ2臨界值，則接受無效假設H0，反之，拒絕無效假設H0。理論χ2臨界值可通過查卡方檢驗臨界值表所得。

3 定軌殘差特征分析

從以上偽距殘差與載波相位殘差中可以看出，CUT0-CO1偽距殘差波形較為雜亂，含有大量隨機噪聲，但仔細分析可以看出較模糊的周期調制信息，CUT0-CO1載波相位殘差波形雜亂，含有大量隨機噪聲，未直觀地觀測出明顯特征。GSMD-CO3偽碼定軌殘差中還有大量隨機噪聲，但波形為均勻分布，GSMD-CO3載波相位定軌殘差較雜亂，但能看見三個明顯的下降峰值。因此，可以看出在北斗定軌殘差中，包含了豐富的信息，從波形分布可以看出定軌殘差并沒有呈現出理想情況下的隨機分布，因此可以推斷在北斗的定軌殘差中包含了可以提取的特征信息，如果能提取這些信息，就能對定軌殘差中的非建模系統誤差進行提取與分離，有望提高北斗衛星定軌的內符合精度。

圖3 CUT0-CO1偽距與載波相位定軌殘差Fig.3 Pseudo-range and carrier phase orbit determination residuals of CUT0-CO1

圖4 GSMD-CO3偽距與載波相位定軌殘差Fig.4 Pseudo-range and carrier phase orbit determination residuals of GSMD-CO3

4 實例分析

下面通過一個CUT0-C01的載波相位定軌殘差數據為例，分析提取定軌殘差中特征信息，然后給出不同接收站針對不同“北斗”GEO衛星定軌殘差中的特征提取結果。

圖5 CUT0-C01載波相位定軌殘差的頻譜分析Fig.5 FFT spectrum of carrier phase orbit determination residuals of CUT0-C01

如第3節所述，CUT0接收站接收C01“北斗”導航衛星的載波相位定軌殘差如圖3所示，對其進行最常用的傅里葉頻譜分析，其頻譜如圖5所示。在圖5中，發現了兩個主要的峰值頻率，分別為1.543×10-5Hz(顯示時，四舍五入到小數點后3位)與3.472×10-5Hz，對應的周期分別為64 809 s與28 802 s?？紤]定軌殘差的物理意義及GEO衛星運動規律，并不能直接得出這兩個頻率的物理含義。由于傳統頻譜分析是基于傅里葉變換，而傅里葉變換是建立在數據平穩性假設條件下的一種頻域全局性變換，對于分析平穩與準平穩時間序列十分有效，但對分析非平穩時間序列的能力有限，不能很好揭示非平穩時間序列的特征[4]，因此需進一步進行“北斗”定軌殘差數據分析。

運用本文提出的總體平均經驗模式分解與Hilbert包絡譜方法對CUT0-C01載波相位殘差進行分析。濾波通帶設置為[3fs/8，fs/2]，濾波器選取為FIR濾波器，期望的信號分解相對誤差最小值e設置為0.002，因此得到EEMD分解方法中的加入白噪聲大小系數為0.026 4，總體平均的次數為174。其中得到8個IMF，和一個殘余項。圖6、圖7是對定軌殘差序列進行EEMD分解后所得的IMF。其中，“signal”表示原始載波相位殘差，“C1”、“C2”、…、“C8”分別為第1至第8個IMF，“r”表示剩余項。可以看出，原始載波相位殘差按照頻帶高低被自動分解為8個IMF和1個剩余項，在此過程中，原始載波相位殘差中的特征信息分別落入8個IMF中。為進一步精確提取這些特征信息，對所得IMF進行譜分析。

圖6 CUT0-C01載波相位定軌殘差EEMD分解1Fig.6 1th EEMD results of carrier phase orbit determination residuals of CUT0-C01

圖7 CUT0-C01載波相位定軌殘差EEMD分解2Fig.7 2nd EEMD results of carrier phase orbit determination residuals of CUT0-C01

將EEMD分解后所得的IMF進行Hilbert包絡譜分析，其結果分別如圖8、圖9所示，分別對應第1個IMF至第8個IMF的Hilbert包絡譜。從圖8中明顯可以找出峰值頻率2.315×10-5Hz；6.944×10-5Hz，從圖9中明顯可以找出峰值頻率1.157×10-5Hz，2.315×10-5Hz，3.472×10-5Hz。

圖8 CUT0-C01定軌殘差EEMD分解后的Hilbert譜1Fig.8 1 th Hilbert spectrum of orbit determination of EEMD results in CUT0-C01

圖9 CUT0-C01定軌殘差EEMD分解后的Hilbert譜2Fig.9 2nd Hilbert spectrum of orbit determination of EEMD results in CUT0-C01

表1 特征頻率與周期對照表

在圖8中1.517×10-5Hz在四舍五入前的真實頻率為1.517 407 407 5×10-5Hz。將特征頻率轉換為時間周期如表1所示。從表1中可以看出，在Hilbert包絡譜中的峰值頻率1.517×10-5Hz嚴格對應86 400 s，正好是1天的積秒。其余峰值特征頻率是1.157×10-5Hz的2,3,6倍頻。因此，初步判定在CUT0-C01的載波相位定軌殘差中存在一個與1天相關的時間周期項，結合“北斗”GEO衛星的軌道運行規律，推斷定軌殘差中可能包含與軌道運動周期嚴格對應的特征信息。

為進一步驗證以上推論，選取CUT0/GSMD/JFNG接收站對多顆“北斗”GEO衛星的定軌殘差運用本文提出的特征提取方法進行分析，找出其中最明顯的兩個特征頻率，其統計結果如表2所示。C01、C03、C04、C05均為北斗GEO衛星。表2中，“/”表示在Hilbert包絡譜中僅存在1個最明顯頻率成分，其他頻率成分相對較弱。從表2中可以看出，在CUT0、GSMD、JFNG站針對多顆北斗GEO衛星的定軌殘差中，均存在與軌道運動周期1天相對應的時間項，表現為存在基頻1.157×10-5Hz以及1.157×10-5Hz的倍頻。初步判斷產生此現象可能的原因為“北斗”導航衛星運動受到太陽光壓的影響產生了周期成分，抑或是“北斗”導航衛星信號從發射到接收過程中經過了較為復雜的傳播路徑，如多徑效應[8]等影響，在接收機端接收到的信號相互調制，產生非線性耦合現象，非線性的主要表現為不同信號頻率成分間的調制現象，產生了頻率邊帶結構。從以上分析可看出，基于總體平均經驗模式分解與Hilbert譜的定軌殘差特征提取方法比傳統傅里葉變換更能準確提取淹沒在噪聲中的特征信息。

表2 定軌殘差中最明顯的頻率特征

圖10 CUT0-C01偽距定軌殘差分解為模型信號與剩余項信號Fig.10 Model signal and residuals signal of Pseudo-range orbit determination residuals of CUT0-CO1

為進一步說明總體平均經驗模式分解方法在處理定軌殘差數據中的有效性，利用卡方檢驗分析處理前后北斗定軌數據的分布特性。以CUT0-C01的偽距定軌殘差數據的EEMD分析為例，通過EEMD分解與重構后，可將CUT0-C01信號分解為一個模型信號與剩余項信號，如圖10所示，模型信號由EEMD分解中的中低頻信號組成，剩余項信號由EEMD分解中的高頻信號組成。

分別對CUT0-C01定軌殘差原始信號與剩余項信號進行卡方檢驗正態分布分析。卡方檢驗的統計自由度為7，顯著性水平為0.05。原始信號所得χ2=157.616 5，剩余項信號所得χ2=5.505 6。經查卡方檢驗臨界值表，對應自由度與顯著性水平的χ2臨界值為14.067 1。因此，運用卡方檢驗方法不難得出，原始信號H0=1，即原始信號不滿足正態分布；剩余項信號H0=0，即剩余項信號滿足正態分布。剩余信號成正態分布，滿足理想定軌殘差的分布特性，推斷定軌殘差中的有效信息主要存在于模型信號中，因此通過卡方檢驗進一步驗證了EEMD方法可有效分離出定軌殘差中包含的特征信息。

5 結束語

1)本文提出了總體平均經驗模式分解與Hilbert譜相結合的“北斗”定軌殘差特征提取方法，準確提取“北斗”導航GEO衛星的非建模系統誤差的特征信息，相比傳統傅里葉頻譜分析方法，在準確提取“北斗”定軌殘差特征信息中具有更高的可靠性和靈敏度。

2)在“北斗”導航衛星的偽距與載波相位的定軌殘差中，存在與軌道運動天周期(86 400 s)相關的特征頻率1.157×10-5Hz及其倍頻。

3)卡方檢驗驗證了總體平均經驗模式分解方法在定軌殘差分析中的有效性。

4)后續需進一步開展的工作包括，一是深入研究“北斗”導航衛星定軌殘差中產生不同周期項的原因，二是對提取出的特征頻率定量信息在定軌殘差中重新建模，去除此項非建模系統誤差，并通過定軌試驗來驗證建模的準確性與可靠性。

致謝

感謝CUT0、GSMD、JFNG等GNSS衛星國際接收站提供的北斗導航衛星的偽距與載波相位觀測量數據。

[1] 郭向, 張強, 趙齊樂, 等.基于單頻星載GPS數據的低軌衛星精密定軌[J]. 中國空間科學技術，2013,33(2):41-46.

GUO XIANG, ZHANG QIANG, ZHAO QILE，et al. Precise orbit determination for LEO satellites using single-frequency GPS observations [J].Chinese Space Science and Technology, 2013, 33(2):41-46.

[2] 施闖, 趙齊樂, 李敏, 等.北斗衛星導航系統的精密定軌與定位研究[J]. 中國科學：地球科學, 2012, 42(6): 854-861.

SHI CHUANG,ZHAO QILE, LI MIN, et al. Precise orbit determination of Beidou satellites with precise positioning [J]. Science China Earth Science., 2012, 42(6)：854-861.

[3] HUANG NORDEN E, ZHENG SHEN, STEVEN R L. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis [J]. Proc. R. Soc. Lond. A, 1998, 454：903-995.

[4] 陳略, 訾艷陽, 何正嘉, 等.總體平均經驗模式分解與1.5維譜的方法研究[J].西安交通大學學報, 2009, 43(5):94-98.

CHEN LUE, ZI YANYANG, HE ZHENGJIA, et al. Research and application of ensemble empirical mode decomposition principle and 1.5 dimension spectrum method [J]. Journal of Xi′an Jiaotong University, 2009, 43(5): 94-98.

[5] CHEN LUE, ZI YANYANG, HE ZHENGJIA, et al. Rotating machinery fault detection based on improved ensemble empirical mode decomposition [J]. Advances in Adaptive Data Analysis, 2014, 6(2):1-24.

[6] WU ZHAOHUA, HUANG NORDEN E. Ensemble empirical mode decomposition: a noise-assisted data analysis method [J]. Advances in Adaptive Data Analysis, 2009, 1(1):1-41.

[7] 何正嘉,訾艷陽,張西寧.現代信號處理技術及應用[M].西安：西安交通大學出版社，2007:61-62.

[8] LIU HUICUI, LI XIAOJING, GE LINLIN, et al. Variable length LMS adaptive filter for pseudorange multipath mitigation [J]. GPS Solutions, 2011, 15(1): 29-38.

(編輯：車曉玲)

Feature Extraction and Distribution Test for Orbit Residuals in Beidou Navigation Satellite

CHEN Lue1,2TANG Geshi1,2CUI Hongzheng1,2CHEN Ming1,2LIU Huicui1,2WANG Mei1,2

(1 Beijing Aerospace Control Center, Beijing 100094)

(2 National Key Laboratory of Science and Technology on Aerospace Flight Dynamic, Beijing 100094)

To analyze the no-modeling system error of orbit residuals in COMPASS, the orbit residuals feature extraction method with the ensemble empirical mode decomposition (EEMD) and Hilbert spectrum was proposed, and the normal distribution character of the orbit residuals chi-square test was analyzed. Firstly, the feature extraction principle of EEMD and Hilbert spectrum was introduced,the EEMD method assisted by the filter was proposed, and the feature extraction model of orbit residuals was established. Secondly, the feature of Pseudo-range and carrier phase orbit residuals in COMPASS was analyzed. Finally, the EEMD and Hilbert spectrum method was applied to analyze the Beidou GEO satellite orbit residuals. The results show that EEMD and Hilbert method can accurately extract the features of 1 day period, which is related to Beidou orbit running period. The chi-square value of the post orbit residuals obtained by EEMD is 5.5, which is smaller than the critical value of chi-square test. It proves that the post orbit residuals are normal distribution. The proposed method is important for separating the no-modeling system error in Beidou navigation satellite orbit residuals to improve the orbit determination precision.

Orbit residuals；Feature extraction；Distribution test；Pseudo-range；Carrier phase；Ensemble empirical mode decomposition；COMPASS

國家自然科學基金(41304026)資助項目

2014-07-11。收修改稿日期：2014-11-06

10.3780/j.issn.1000-758X.2015.02.001

陳 略 1983年生，2009年獲西安交通大學機械工程及自動化專業碩士學位，工程師。研究方向為無線電測量與科學應用、地球自轉參數解算與預報。