一種基于小波分析的衛(wèi)星鐘差數(shù)據(jù)粗差處理方法

王 威 許 芬 王宇譜

1 北京衛(wèi)星導(dǎo)航中心，北京市北清路22號，100094 2 北方工業(yè)大學(xué)電氣與控制工程學(xué)院，北京市晉元莊路5號，100144

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

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

1 算法原理

小波分析理論具有較好的局部特征分辨能力，本文利用小波分析的優(yōu)點，結(jié)合衛(wèi)星鐘差數(shù)據(jù)及其異常點的特性，給出一種基于小波分析的衛(wèi)星鐘差數(shù)據(jù)預(yù)處理方法。

1.1 基于小波分析的鐘差數(shù)據(jù)預(yù)處理原理

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

本文基于小波分析設(shè)計的衛(wèi)星鐘差數(shù)據(jù)粗差處理原理與過程見圖1。

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

圖1 數(shù)據(jù)處理流程Fig.1 Data process procedure

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

1.2 小波函數(shù)與分解尺度的選型

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

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

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

1.3 不同分解尺度處理鐘差數(shù)據(jù)的差異分析

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

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

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

圖3和4分別為C09衛(wèi)星鐘差頻率數(shù)據(jù)經(jīng)過db3小波用分解尺度2和3進(jìn)行分解后各層小波系數(shù)及粗差剔除情況?？梢郧逦乜吹?，數(shù)據(jù)中存在粗差值，且在各層高頻小波系數(shù)中能明顯觀察到系數(shù)發(fā)生跳變。用尺度2分解后，在低頻小波系數(shù)中也能明顯觀察到粗差點；而用尺度3分解后，低頻小波系數(shù)無明顯的粗差點，與§1.2關(guān)于分解尺度越大低頻小波系數(shù)分解得到的信息越少的分析結(jié)論相互印證。因此，選擇合適尺度進(jìn)行分解對粗差剔除的效果有直接的影響。經(jīng)過處理，各層小波系數(shù)的跳變值被剔除，通過小波重構(gòu)則可得到剔除粗差后的鐘差頻率數(shù)據(jù)。

1.4 兩種數(shù)據(jù)預(yù)處理方法的差異分析

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

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

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

通過對小波函數(shù)與分解尺度的選型進(jìn)行分析發(fā)現(xiàn)，選用分解尺度3或4的db3小波函數(shù)進(jìn)行預(yù)處理能獲得較好的結(jié)果。因此，本文在隨后的實驗中統(tǒng)一使用分解尺度3的db3小波函數(shù)。

2 實驗分析

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

2.1 數(shù)據(jù)擬合情況對比

利用小波分析對CODE解算的BDS各顆衛(wèi)星的鐘差數(shù)據(jù)進(jìn)行預(yù)處理，分別比較各顆衛(wèi)星數(shù)據(jù)跨度為24 h、48 h和72 h的擬合精度。

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

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

2.2 數(shù)據(jù)預(yù)報情況對比

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

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

表1 2種預(yù)處理方法在不同數(shù)據(jù)長度時的預(yù)報精度對比

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

3 結(jié) 語

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