一種基于小波分析的衛星鐘差數據粗差處理方法

王 威 許 芬 王宇譜

1 北京衛星導航中心，北京市北清路22號，100094 2 北方工業大學電氣與控制工程學院，北京市晉元莊路5號，100144

衛星鐘差數據能直接反映星載原子鐘的工作狀態，作為衛星導航系統不可或缺的一種數據產品，其質量的優劣對導航系統的性能至關重要。由于在軌衛星容易受外部環境變化等因素的影響，導致獲取的衛星鐘差數據經常出現間斷、跳變和粗差等異常情況[1]，而含有異常數據的衛星鐘差不僅會降低衛星鐘性能分析的可靠性，還會導致鐘差預報的精度隨預報時長的增加快速衰減。由此可見，衛星鐘差數據的預處理是原子鐘穩定性分析、鐘差預報等科學研究和工程應用的前提和基礎。針對鐘差數據的粗差探測與剔除等問題，較為常用的方法是中位數法(MAD)及其改進方法[1-5]。但MAD法參數的選取是基于經驗值，若參數選取較大則粗差剔除得不夠徹底；若參數選取過小，則部分有效數據會被當作粗差剔除。另外，還有其他的一些預處理方法[6-10]。

基于此，本文提出一種基于小波分析的數據預處理方法，能充分利用小波分析對數據局部特征的分辨能力，有效處理小量級的粗差值，同時顧及預處理后對鐘差預報精度的影響，保證經過預處理后的數據能獲得較高精度的預報參數。使用CODE解算的BDS衛星鐘差數據進行數據預處理及擬合、預報，并與常用的MAD法進行對比分析，結果表明，新方法能夠有效剔除數據中的異常值，得到較好的預處理數據，且基于新方法處理后的鐘差數據能夠有效提高衛星鐘差預報的精度。

1 算法原理

小波分析理論具有較好的局部特征分辨能力，本文利用小波分析的優點，結合衛星鐘差數據及其異常點的特性，給出一種基于小波分析的衛星鐘差數據預處理方法。

1.1 基于小波分析的鐘差數據預處理原理

利用小波分析進行粗差探測的主要依據在于數據序列中粗差點相較鄰近數據有較大起伏，并且這些起伏的特征在不同尺度變換的小波系數中均有對應。因此，可以對數據進行小波分解，并通過處理小波系數的起伏情況進行粗差探測和剔除。

本文基于小波分析設計的衛星鐘差數據粗差處理原理與過程見圖1。

首先對鐘差頻率數據進行小波分解，得到低頻小波系數cai和各層高頻小波系數cdj，再對分解所得的低頻、高頻小波系數進行粗差處理。需要注意的是，在粗差處理過程中，對處理效果影響最大的因素是低頻小波系數cai，因此采用低頻小波的中位數作為閾值，即med=median{|cai|}，若cai>med/3，則將該點的小波系數設為粗差點[10]。對于高頻小波系數cdj，利用σ=median{|cdj|/0.674 5}計算各層小波系數的方差，將|cdj|>3σ的小波系數點判定為粗差點。考慮到個別衛星的頻率數據較為明顯地偏離零值，若直接將粗差值置零，則置

圖1 數據處理流程Fig.1 Data process procedure

零點相當于成為了新的粗差值；若簡單地將粗差值剔除，則會引起不必要的數據缺失。因此，本文使用三次分段樣條方法對探測出的粗差點進行內插并替換，以減小對鄰近數據的影響。將處理后的小波系數進行重構，獲得粗差處理后的鐘差頻率數據，同時利用參數為5的MAD法對原始頻率數據進行剔除粗差處理，分別將2種粗差預處理后的頻率數據恢復為相位數據，并進行擬合與預報，對比分析二者的差異。

1.2 小波函數與分解尺度的選型

利用小波函數進行鐘差數據預處理時，可選用多種小波函數，常用的有db系、sym系和coif系，每種小波函數又可搭配多種分解尺度。因此，衛星鐘差數據預處理的效果同時受小波函數和分解尺度兩個因素的影響，選取最適合處理衛星鐘差數據的小波函數和分解尺度對數據的擬合預報至關重要。文獻[3]和[10]分析了不同小波函數處理BDS鐘差數據的差異，認為db系小波的處理精度要高于sym系和coif系小波，因此，本文選用分解尺度為2、3、4、5、6的dbN(N=2,3,4,5,6,7,8,9,10)小波分別處理衛星鐘差。考慮到C10衛星的數據相對完整連續，此處對該衛星的數據進行預處理(選用其他衛星亦可)，并利用二次多項式模型對數據進行擬合，統計擬合殘差的RMS值，隨后利用二次多項式模型對數據進行2 h預報，統計預報殘差的RMS值。

圖2中縱軸表示不同小波函數的擬合預報精度，橫軸表示不同分解尺度的擬合預報精度，其中藍色越深表示擬合或預報精度越高，紅色越深表示精度越低。對于同一分解尺度而言，不同小波處理的擬合精度相當(除分解尺度5的db10小波函數擬合精度較差及分解尺度6的db5～10小波外)，大部分小波函數經過分解尺度6的小波分解處理后擬合精度不高，其原因是當分解尺度越大時，低頻小波系數分解得到的信息越少，導致重構后的數據存在失真。對于同一小波函數而言，隨著分解尺度的增加，擬合殘差整體呈逐漸變小的趨勢，其原因是當分解尺度增大后，對數據去噪和平滑的效果更為明顯，處理后數據的高頻噪聲被顯著消除，離散程度變小，其擬合殘差也相應變小。擬合精度最高的組合是分解尺度6的db3小波函數。對于預報而言，分解尺度為2、3、4的各小波均有較高的預報精度，其中分解尺度3的小波預報精度整體較高。另外，在分解尺度大于4的情況下，預報精度反而出現衰減。從同一小波函數來看，各小波函數在分解尺度小于4時，均有較好的預報精度。從整體來看，db2、db3、db4小波的預報結果較好。

圖2 不同小波函數和分解尺度各層小波系數及粗差剔除情況Fig.2 Wavelet coefficients and gross error elimination at different wavelet and decomposition scales

1.3 不同分解尺度處理鐘差數據的差異分析

基于小波分析的數據預處理受所選小波函數及分解尺度的影響，其中分解尺度對預處理的影響更為顯著，因此應著重對分解尺度的影響展開分析。

圖3 尺度2的小波分解后的各層小波系數及粗差剔除情況Fig.3 Wavelet coefficients and gross error elimination of each layer after 2td dimension wavelet decomposition

圖4 尺度3的小波分解后的各層小波系數及粗差剔除情況Fig.4 Wavelet coefficients and gross error elimination of each layer after 3th dimension wavelet decomposition

圖3和4分別為C09衛星鐘差頻率數據經過db3小波用分解尺度2和3進行分解后各層小波系數及粗差剔除情況。可以清晰地看到，數據中存在粗差值，且在各層高頻小波系數中能明顯觀察到系數發生跳變。用尺度2分解后，在低頻小波系數中也能明顯觀察到粗差點；而用尺度3分解后，低頻小波系數無明顯的粗差點，與§1.2關于分解尺度越大低頻小波系數分解得到的信息越少的分析結論相互印證。因此，選擇合適尺度進行分解對粗差剔除的效果有直接的影響。經過處理，各層小波系數的跳變值被剔除，通過小波重構則可得到剔除粗差后的鐘差頻率數據。

1.4 兩種數據預處理方法的差異分析

為比較本文提出的數據預處理方法和MAD法在粗差探測方面的差異，選用分解尺度3的db3小波函數作為小波分析的數據預處理方法，并與參數為5的MAD法進行對比，通過人工判讀的方式評估粗差探測的效果。

圖5(a)反映的是包含粗差的原始鐘差頻率，圖5(b)和5(c)分別為通過MAD法和小波分析預處理后的頻率數據。MAD法和小波分析均能有效地剔除粗差數據，說明小波分析處理粗差值具有合理性。對比MAD法和分解尺度3小波函數的處理結果可以看到，經參數為5的MAD法處理后，較大的粗差值被有效探測并剔除，但仍存有一部分較小的頻率跳變值未能被探測到；而使用小波分析處理后，部分高頻噪聲被消除，量級較小的粗差值會隨著高頻噪聲的消除被一并剔除。從粗差探測與剔除效果來看，本文方法處理后的數據比MAD法更為集中和平滑。選用分解尺度較大的小波則集中程度更為明顯，這在一定程度上改變了原始鐘差數據，對于評估衛星鐘的性能而言，勢必會降低評估結果的客觀性，但對于鐘差預報而言，則有利于參數的估計。由于衛星鐘差數據在導航定位的應用中是需要進行預報的，因此預報精度的提升對于導航定位有著十分實用的價值。據此，本文設計了一套基于鐘差預報的量化實驗，將MAD法和小波分析處理后的頻率數據恢復為相位數據，并對該數據進行擬合與預報，對比分析2種預處理方法的預報精度。

圖5 MAD法與小波分析預處理后的對比Fig.5 Comparison between MAD method and wavelet analysis

通過對小波函數與分解尺度的選型進行分析發現，選用分解尺度3或4的db3小波函數進行預處理能獲得較好的結果。因此，本文在隨后的實驗中統一使用分解尺度3的db3小波函數。

2 實驗分析

本文使用2018-02-01～02-10共10 d的CODE解算的BDS衛星鐘差相位數據(ftp://cddis.gsfc.nasa.gov/pub/gps/products/mgex/)，數據采樣間隔為300 s。

2.1 數據擬合情況對比

利用小波分析對CODE解算的BDS各顆衛星的鐘差數據進行預處理，分別比較各顆衛星數據跨度為24 h、48 h和72 h的擬合精度。

由圖6可知，尺度3的db3小波分解和MAD法預處理后數據的擬合精度相當，并且擬合殘差隨著數據跨度的增加而變大。分析其原因，一方面是由于多星定軌解算的鐘差存在跨天跳變，導致跨度較長的數據存在不連續的情況，影響擬合殘差的統計結果；另一方面，跨度較長的數據包含變頻和跳相的可能性會增加，最終影響擬合精度。

圖6 db3小波處理不同數據跨度的擬合情況Fig.6 Fitting residuals of db3 wavelet preprocess with different data spans

2.2 數據預報情況對比

分別統計使用小波分析和MAD法預處理的擬合資料跨度為24 h、48 h和72 h的數據2 h的預報精度，并比對2種方法以24 h擬合資料分別預報2 h、4 h、6 h和12 h的預報精度。為了更集中、綜合地評價MAD法和基于小波分析的數據預處理方法在鐘差數據預報方面的效果，對計算結果進行算術平均，得到的預報精度不是比較某顆衛星的預報效果，而是比較2種預處理方法的效果。

表1統計了2種預處理方法在擬合資料長度不同的情況下2 h的預報精度對比，可以看出，2種預處理方法2 h的預報精度均優于2 ns，基于小波分析的預處理方法得到的數據在鐘差預報方面要優于MAD法，預報精度至少提高10%。表2統計了2種預處理方法在擬合資料長度相同的情況下不同時長的預報精度的對比，可以看出，預報精度隨預報時長的增加而增大，基于小波分析預處理方法的預報精度比MAD法有提升，但提升幅度隨預報時長的增加而降低。

表1 2種預處理方法在不同數據長度時的預報精度對比

表2 2種預處理方法的不同時長的預報精度對比

3 結 語

本文提出一種基于小波分析的衛星鐘差數據預處理方法，通過合理設置各層小波系數的異常數據剔除閾值，可有效剔除數據中的粗差點。分析使用不同分解尺度、不同小波函數對粗差處理的影響后發現，分解尺度3和4的db3小波函數最適合進行粗差處理。另外，本文還進行了不同數據跨度的擬合精度及預報精度分析，從結果上看，擬合精度隨分解尺度的增大而逐漸提高。對比MAD法，小波分析的預處理方法在2 h短期預報方面有較大優勢，精度提升超過10%。