梁益豐，許江寧，吳 苗，何泓洋

（海軍工程大學 電氣工程學院，武漢 430033）

全球衛(wèi)星導航系統(tǒng)（Global Navigation Satellite System,GNSS）本身是高精度時間同步系統(tǒng)，實現(xiàn)定位與授時需要穩(wěn)定的時間基準和高精度鐘差預報技術(shù)。我國北斗衛(wèi)星導航系統(tǒng)（BeiDou Navigation Satellite System,BDS）于2020年7月全面開通，隨著衛(wèi)星數(shù)量和性能的不斷提升，將衛(wèi)星原子鐘逐步納入系統(tǒng)時間基準計算，最終實現(xiàn)完全基于衛(wèi)星原子鐘的系統(tǒng)時間基準生成，對于BDS安全性和性能提升有重要意義[1]，因此有必要開展基于衛(wèi)星鐘差數(shù)據(jù)的時間尺度算法、鐘差預報等相關(guān)技術(shù)研究。

由于運行環(huán)境復雜，衛(wèi)星鐘差數(shù)據(jù)通常包含趨勢分量、多個周期分量、隨機分量等成分。常用的最小二乘擬合趨勢項方法不夠準確，往往導致頻率計算偏差，影響鐘差預報精度和ALGOS類時間尺度算法中的頻率預測[2]；對擬合殘差進行頻譜分析時，由于噪聲項和殘留趨勢項的干擾，可能導致周期辨識不準、幅值產(chǎn)生偏差[3]；利用鐘差或頻率數(shù)據(jù)計算穩(wěn)定度時，周期項通常會造成Allan偏差出現(xiàn)異常波動，不利于穩(wěn)定度分析與噪聲系數(shù)擬合，進而影響時間尺度算法中的權(quán)重選取[4]。因此，部分信號處理技術(shù)被先后應(yīng)用于鐘差數(shù)據(jù)分析：傅里葉變換將信號從時域轉(zhuǎn)換到頻域，以分離和區(qū)分平穩(wěn)信號和噪聲，但對于非平穩(wěn)和非線性信號無效；小波變換可以在時域和頻域中可視化信號，并且對非平穩(wěn)信號有效，但其假設(shè)信號在小波窗口中平穩(wěn)，且存在小波基選擇問題[5]；Kalman濾波被廣泛用于鐘差去噪和狀態(tài)估計，但難以同時分解各種成分，且建模和實現(xiàn)較為復雜[6]。

經(jīng)驗模態(tài)分解（Empirical Mode Decomposition,EMD）是黃鍔博士提出的一種自適應(yīng)信號時頻處理方法，適用于非線性非平穩(wěn)信號的辨識、去噪與預測[7]，該方法及其衍生模型近些年逐漸被應(yīng)用于鐘差數(shù)據(jù)處理。在頻率穩(wěn)定度方面，朱江淼利用集合經(jīng)驗模態(tài)分解（Ensemble Empirical Mode Decomposition,EEMD）對原子鐘頻差數(shù)據(jù)進行降噪，效果優(yōu)于小波閾值方法[8]；惠恬將EMD與小波降噪相結(jié)合，一定程度上提升了頻率穩(wěn)定度[9]；Aly I.Mostafa利用EMD提取隨機分量，分別提升了單鐘和鐘組時間尺度的短期穩(wěn)定度[10]。在周期項分析方面，伍貽威提出了原子鐘模型和頻率穩(wěn)定度的系統(tǒng)性分析方法，并應(yīng)用于地面氫鐘數(shù)據(jù)分析[11]；肖勝紅提出了一種奇異譜分解與傅里葉帶通濾波器相結(jié)合的周期項提取方法，提高了24 h周期項提取效率[12]；李驍逸采取頻譜分析確定周期大小與幅值，對提取周期項的鐘差數(shù)據(jù)進行穩(wěn)定度分析，表明周期項對頻率的短期穩(wěn)定度有顯著影響[1]。上述分析大多集中于對鐘差某一種成分的判定，或通過較多運算逐步分析各成分特征。同時，常用EMD方法存在模態(tài)混疊問題，在EMD基礎(chǔ)上加入成對正負高斯白噪聲的EEMD可以削弱模態(tài)混疊現(xiàn)象，然而其分量總會殘留一定白噪聲，影響后續(xù)分析與處理[13]。為此，TORRES等提出自適應(yīng)噪聲完備集合經(jīng)驗模態(tài)分解（Complete Ensemble Empirical Mode Decomposition with Adaptive Noise,CEEMDAN）方法，在每次求解分量后重新給殘值加入白噪聲，然后逐次迭代求解[14]，較好改善了模態(tài)混疊現(xiàn)象。

針對上述問題，本文提出一種能夠同步提取鐘差不同成分的“分解—辨識—重構(gòu)”方法，在分析鐘差主要特征的基礎(chǔ)上，使用CEEMDAN方法完成數(shù)據(jù)分解，綜合應(yīng)用排列熵算法與t檢驗分析原理，對信號分量進行成分辨識，進而重構(gòu)得到趨勢分量、周期分量、隨機分量，通過BDS鐘差實測數(shù)據(jù)驗證了方法的有效性和實用性。

1 BDS衛(wèi)星原子鐘鐘差特性分析

1.1 衛(wèi)星原子鐘時差模型

以相位、頻率、頻漂三種參數(shù)組成的多項式模型是鐘差的主要成分，也通常被應(yīng)用于擬合趨勢分量，周期分量與隨機分量與趨勢項相比數(shù)量級較小，但對穩(wěn)定度的影響會逐漸積累。

1.2 各分量對穩(wěn)定度的影響

以Allan方差為例，分析各分量在頻率穩(wěn)定度中的影響。文獻[15]詳細推導了噪聲擴散系數(shù)與Allan方差的關(guān)系：

其中右邊第1項為觀測噪聲，第4項為頻漂，二者在對數(shù)Allan方差圖中的斜率分別為-2和2。當周期表征較為明顯時，Allan方差也將出現(xiàn)周期波動，表現(xiàn)為在某一平滑時間內(nèi)的異常凸起。根據(jù)以上分析，常規(guī)噪聲和頻漂對頻率穩(wěn)定度的影響易于判斷，然而衛(wèi)星原子鐘的多周期項會導致頻率穩(wěn)定度出現(xiàn)不確定因素，進而影響擴散系數(shù)擬合和性能分析。顯然，將衛(wèi)星鐘差進行準確分解與重構(gòu)十分必要，由于攜帶周期性波動的衛(wèi)星鐘差進行系統(tǒng)時間基準計算會將上述波動引入系統(tǒng)時間，因此未消除周期項的鐘差數(shù)據(jù)不應(yīng)納入系統(tǒng)時間基準計算。

2 基于CEEMDAN分解的信號分析方法

2.1 CEEMDAN分解原理

計算去除第二個模態(tài)分量后的殘差：

重復上述步驟，直到獲得的殘差信號為單調(diào)函數(shù)，不能繼續(xù)分解，算法結(jié)束。此時得到的本征模態(tài)分量數(shù)量為K，則原始信號()x T被分解為：

由于加入白噪聲方式的改進，CEEMDAN在較小的平均次數(shù)下就能達到極小的重構(gòu)誤差，完備性優(yōu)于EEMD方法，也因此具備更少的平均次數(shù)和更快的計算速度[14]。此外，EEMD分解還可能出現(xiàn)多個幅值很小的低頻IMF分量，CEEMDAN方法能夠改善此現(xiàn)象。其具體分解流程如圖1所示。

圖1 CEEMDAN分解流程Fig.1 The CEEMDAN decomposition process

2.2 排列熵算法

排列熵算法可以衡量系統(tǒng)復雜度，檢測時間序列隨機項突變的方法，計算簡單快速，魯棒性強，已被廣泛應(yīng)用于非線性數(shù)據(jù)處理與分析[16]。第i個模態(tài)分量IMFi(T)的排列熵值求解過程如下：

步驟1：對長度為N的IMFi(T)進行m維相空間重構(gòu)，生成K×m的矩陣IMF：

排列熵值大小表示時間序列隨機程度：熵值越小，說明時間序列越簡單、規(guī)則；熵值越大，則時間序列越復雜、隨機。由此，引用該算法定性分析鐘差信號中的信號主導分量和噪聲主導分量。

2.3 模態(tài)分量特征分析

盡管排列熵算法能夠有效辨識信號與噪聲主導分量，但因為衛(wèi)星工作環(huán)境復雜，經(jīng)分解的鐘差數(shù)據(jù)臨界模態(tài)仍然可能存在多種信號。為此，引用文獻[17]所采取的t檢驗方法對模態(tài)分量進行特征分析，以確保重構(gòu)信號的準確性。

t檢驗是用t分布理論來推論差異發(fā)生的概率，從而比較兩個平均數(shù)的差異是否顯著。分析CEEMDAN原理可知，其IMF分量應(yīng)滿足上、下包絡(luò)線相對于時間軸局部對稱。則高頻IMF分量的上下包絡(luò)線基本由眾多的信號峰值點連接得到，對稱的包絡(luò)線意味著IMF數(shù)據(jù)基本對稱，數(shù)據(jù)均值趨近于0；低頻IMF信號周期大，包絡(luò)線由少量峰值插值獲取，與原信號趨勢走向關(guān)聯(lián)度較低，所以信號分量并不對稱，很難保證均值為0。由此，將分解得到的K個IMF做高低頻區(qū)分。不妨假設(shè)IMF1為指標1，IMF1+IMF2為指標2，以此類推，前i個IMF的和為指標i，自指標1開始計算均值，對其是否顯著區(qū)別于0進行t檢驗分析，當顯著區(qū)別于0時停止計算。t檢驗統(tǒng)計量為[17]：

2.4 鐘差數(shù)據(jù)分解與重構(gòu)流程

綜合衛(wèi)星鐘差數(shù)據(jù)處理、信號分解、模態(tài)分量特征分析等過程，給出鐘差數(shù)據(jù)分解與重構(gòu)流程如圖2所示，具體步驟如下：

圖2 鐘差數(shù)據(jù)分解與重構(gòu)流程圖Fig.2 Flow chart of clock difference data decomposition and reconstruction

步驟1：鐘差數(shù)據(jù)預處理，利用中位數(shù)探測法完成粗差剔除，并進行平滑插補。考慮到頻率數(shù)據(jù)有效位次多于鐘差數(shù)據(jù)，因此做一次差分求平均頻率()y T；

步驟2：將()y T通過CEEMDAN分解為K個本征模態(tài)分量，按照高頻到低頻順序排列；不失一般性，分解時設(shè)置附加噪聲標準差與()y T標準差之比為0.2，信號平均次數(shù)為100次；

步驟3：計算各個模態(tài)的排列熵值，根據(jù)相鄰值的變化確定信號主導分量和噪聲主導分量（一般第K個模態(tài)為趨勢項，此處通過排列熵進行核驗）；

步驟4：對前M(M≤K)模態(tài)之和進行t檢驗，根據(jù)其均值是否顯著區(qū)別于0判斷高頻和低頻分量；

步驟5：當排列熵確定的噪聲主導分量和t檢驗得到的高頻分量相同時，直接判定為鐘差隨機分量；若存在辨識不一致分量，則繼續(xù)分解相應(yīng)模態(tài)，直至達到一致，以確保周期信號與噪聲信號完全分離；

步驟6：對被判定為隨機分量的模態(tài)累加得到隨機項，除隨機項和趨勢項之外的模態(tài)累加得到周期項，數(shù)據(jù)分解與重構(gòu)完成。

3 算例分析

采用德國地學研究中心（Deutsches geoforschungs zentrum,GFZ）發(fā)布的BDS衛(wèi)星精密鐘差數(shù)據(jù)，選擇數(shù)據(jù)連續(xù)性較好的C04、C13、C37、C40衛(wèi)星原子鐘，包含GEO（C04）、MEO（C37）、IGSO（C13與C40）軌道類型，銣原子鐘（C04、C13、C37）和氫原子鐘（C40）鐘型，取2021年7月4日至7月13日共計10天2880個數(shù)據(jù)點，歷元間隔為5 min。

3.1 鐘差數(shù)據(jù)分解與重構(gòu)效果

以C37為例展示數(shù)據(jù)分解與重構(gòu)的完整過程，其平均頻率數(shù)據(jù)經(jīng)CEEMDAN分解得到11個IMF，各分量與原始數(shù)據(jù)圖如圖3所示。原始信號與IMF之和做差，其殘差均方根為6.24 ×10-28，說明C37平均頻率數(shù)據(jù)得到了完全分解。計算各模態(tài)排列熵值，IMF11計算結(jié)果為0.0820，表明幾乎不含波動分量，因此直接作為趨勢項。其余分量排列熵如圖4所示。

圖3 C37衛(wèi)星原子鐘頻率數(shù)據(jù)CEEMDAN分解結(jié)果Fig.3 CEEMDAN decomposition results for the C37 satellite clock frequency data

圖4 C37衛(wèi)星原子鐘各頻率模態(tài)排列熵值Fig.4 The entropy of each mode for C37

由排列熵顯著拐點初步判斷，第1~6個IMF為噪聲主導分量，第7~10個IMF為信號主導分量。此外，t檢驗也在指標7處顯著不為0，據(jù)此分別將第1~6個IMF相加得到高頻分量、第7~10個IMF相加得到低頻分量，完成信號重構(gòu)如圖5所示。

圖5 C37衛(wèi)星原子鐘各分量重構(gòu)結(jié)果Fig.5 Results of each component reconstruction for C37

在分析過程中，C04原子鐘經(jīng)一次分解后的排列熵與t檢驗在8IMF分量處結(jié)果不同、進行了二次分解，一方面反映出排列熵和t檢驗結(jié)果大多相符，另一方面也說明了本文綜合分析策略的合理性與必要性。4顆衛(wèi)星原子鐘主要成分重構(gòu)結(jié)果如圖6所示。

圖6 4顆衛(wèi)星原子鐘各分量重構(gòu)結(jié)果Fig.6 Results of each component reconstruction for four clocks

通過圖6可以得到，4顆原子鐘頻率變化趨勢分別為遞增（C04）、遞減（C37）與近似平穩(wěn)（C13、C40），CEEMDAN方法能夠準確分解得到各種變化特點的趨勢項，即使近似平穩(wěn)的微小波動也能跟蹤；各原子鐘都含有多個規(guī)律周期，噪聲項呈無規(guī)律變化，數(shù)值普遍大于周期項。不同原子鐘初步對比結(jié)果表明，C37原子鐘頻率值小于C04與C13，但周期分量顯著高于其余衛(wèi)星原子鐘，表明MEO衛(wèi)星受周期項影響更為明顯；C40原子鐘頻率值、周期項、噪聲項均為最小，反映了氫原子鐘的優(yōu)良性能。

3.2 周期項特征分析

現(xiàn)有研究已經(jīng)對GNSS衛(wèi)星原子鐘長期特性進行了較全面的分析，表明其包含與系統(tǒng)軌道周期密切相關(guān)的周期項。例如，BDS衛(wèi)星原子鐘的主周期項通常有24 h、12 h、8 h、6 h，Galileo衛(wèi)星原子鐘的主周期項約為14 h和7 h[18,19]等（軌道周期14.039 h）。當特性分析時間較長時，隨時間累積的周期項通常在頻譜圖中呈現(xiàn)較明顯的波峰，但當特性分析時間較短時，鐘差擬合殘差和噪聲可能會導致頻譜圖出現(xiàn)局部峰值。以上結(jié)論可為頻譜分析提供一定參考。

通過常用二次多項式擬合相位得到相位擬合殘差，以CEEMDAN分解重構(gòu)得到頻率周期數(shù)據(jù)，然后分別進行快速傅里葉變換獲得衛(wèi)星原子鐘的周期特性。4顆衛(wèi)星原子鐘頻譜圖如圖7所示，可以看出，BDS各類衛(wèi)星原子鐘均存在多種周期項，對比擬合殘差和提取周期項頻譜分析結(jié)果，主要結(jié)論如下：

圖7 4顆衛(wèi)星原子鐘周期項分析結(jié)果Fig.7 Analysis of periodic term of four clocks

（1）二次多項式擬合方法難以完整提取鐘差趨勢項、同時包含大量高頻分量，其殘存趨勢分量導致頻譜圖靠近y軸處出現(xiàn)不規(guī)則峰值，高頻分量導致主周期項右側(cè)出現(xiàn)大量不規(guī)則波動；本文方法提取的周期信號，其低頻部分未出現(xiàn)不規(guī)則峰值、主周期項特征明顯，且高頻處基本沒有異常波動，表明提取所得周期項基本不含趨勢和隨機分量，提取效果較好；

（2）主周期變化方面，C04原子鐘擬合殘差與所提取周期項在峰值數(shù)量與主周期方面區(qū)別較大，其擬合殘差周期項排序依次為24 h、12 h、6 h、8 h、4.8 h、4 h、3.43 h等，所提取周期項主周期為12 h、24 h、8 h、6 h（前兩項幅值接近），由于周期主要因衛(wèi)星運行導致、擬合殘差的高頻周期項數(shù)值大多是主周期項的公約數(shù)，合理推測：擬合殘差的高頻周期項可能為周期信號在隨機分量干擾下耦合產(chǎn)生。而本文方法有效分離了周期信號與隨機信號，避免了耦合現(xiàn)象；

（3）不同類型衛(wèi)星原子鐘的主周期因軌道區(qū)別而有所差異，GEO衛(wèi)星鐘主周期為12 h、24 h、8 h、6 h，MEO衛(wèi)星鐘主周期為12 h、8 h、6 h。對于IGSO衛(wèi)星，以擬合殘差分析得到的C13原子鐘主周期為24 h、6 h，C40原子鐘主周期為24 h、8 h，根據(jù)周期項產(chǎn)生原因，相同軌道的衛(wèi)星原子鐘周期項通常一致，說明二者擬合殘差數(shù)據(jù)中的多余信號干擾了頻譜分析結(jié)果，部分周期項被淹沒，這也是6 h、8 h波峰附近出現(xiàn)較多不規(guī)則信號的主要原因。以所提取的周期項分析時，二者主周期項均為24 h、8 h、6 h，其中C13所提取6 h周期項左側(cè)仍然存在多余信號，可能是由多周期項功率譜旁瓣所致，但相比擬合殘差頻譜圖已有明顯改善，提取效果將在3.3節(jié)得到進一步驗證。此外，C37原子鐘周期項幅值占比最大、C40原子鐘周期項幅值占比最小，反映出MEO衛(wèi)星受到周期項干擾最大，該現(xiàn)象可能是由于MEO衛(wèi)星與地球運行非同步所致；GEO和IGSO衛(wèi)星受周期項影響相對較小，其中IGSO氫原子鐘最小。

3.3 頻率穩(wěn)定度分析

結(jié)合衛(wèi)星原子鐘時差模型和各分量特征，對BDS衛(wèi)星原子鐘的頻率穩(wěn)定度展開研究，以進一步檢驗本文方法提取周期項的準確性，并定量分析周期項對衛(wèi)星原子鐘穩(wěn)定度的影響。分別計算4顆衛(wèi)星原子鐘原始數(shù)據(jù)、周期項數(shù)據(jù)、去除周期項數(shù)據(jù)的Allan偏差，如圖8所示。可以看出，4顆衛(wèi)星鐘原始穩(wěn)定度曲線都存在異常凸起，與所提取周期項的穩(wěn)定度曲線吻合，說明本文所提取的周期項能夠準確表征對穩(wěn)定度的影響，在去除周期項后，各原子鐘在不同平滑時間內(nèi)的穩(wěn)定度指標均得到一定提高。對于不同軌道類型的原子鐘，MEO衛(wèi)星原子鐘周期項對穩(wěn)定度的影響最明顯，峰值約為2.21×10-13，IGSO氫原子鐘周期項影響最小，峰值約為3.02 ×10-14；MEO衛(wèi)星鐘周期項對萬秒穩(wěn)的影響最為明顯，其余衛(wèi)星鐘的20000~40000 s穩(wěn)定度受周期項影響最大，與相應(yīng)主周期項顯著相關(guān)。但是不論何種軌道類型、何種原子鐘類型的衛(wèi)星，校正后的穩(wěn)定度曲線都顯著改善了原穩(wěn)定度曲線的異常凸起。

圖8 4顆衛(wèi)星原子鐘周期項分析結(jié)果Fig.8 Analysis of periodic term for four clocks

將周期項提取與穩(wěn)定度分析結(jié)果總結(jié)如表1所示，所有衛(wèi)星的頻率穩(wěn)定度性能都獲得了顯著的提升，萬秒頻率穩(wěn)定度提升幅度依次為14.0%、47.2%、17.6%、7.8%，平均為21.6%。其中，C37衛(wèi)星提升效果最為顯著，萬秒頻率穩(wěn)定度由2.44 ×10-13提升到1.29 ×10-13，其余衛(wèi)星原子鐘穩(wěn)定度提升最顯著的采樣時間主要在20000~40000 s之間，因此萬秒穩(wěn)改善程度相對不明顯。對于MEO衛(wèi)星數(shù)量較多的BDS-3系統(tǒng)，這種改進明顯有利于時間基準的建立與維持。

表1 衛(wèi)星原子鐘數(shù)據(jù)綜合分析結(jié)果Tab.1 Comprehensive analysis of satellite clock data

4 結(jié)論

分析了衛(wèi)星鐘差不同分量對于頻率穩(wěn)定度的影響，闡述了鐘差數(shù)據(jù)分解與重構(gòu)的必要性，提出了綜合CEEMDAN分解、排列熵原理和t檢驗的BDS衛(wèi)星鐘差信號分解與周期項提取方法?；?顆不同類型BDS衛(wèi)星原子鐘的實測數(shù)據(jù)展開分析，結(jié)果表明，本文方法能夠有效、準確地分解衛(wèi)星原子鐘趨勢分量、周期分量與隨機分量，便于各種成分的定量分析和原子鐘性能評估；所提取周期分量的頻譜圖相比常規(guī)多項式擬合殘差頻譜圖更為清晰、周期特征更加明顯；對提取周期分量進行穩(wěn)定度分析，其穩(wěn)定度曲線與鐘差穩(wěn)定度曲線的異常凸起吻合度較高，進一步驗證了本文方法的有效性。去除周期分量后，4顆BDS衛(wèi)星原子鐘的萬秒穩(wěn)平均提升21.6%，將有助于高穩(wěn)定度時間基準的生成與維持。

實測數(shù)據(jù)分析部分主要展示了BDS衛(wèi)星鐘差分解與周期項提取效果，但通過研究過程不難看出，基于鐘差信噪分離的思路與方法在頻率和鐘差預報、原子鐘數(shù)據(jù)降噪、原子時算法等方面都有廣闊的應(yīng)用前景[20]。