高性能原子鐘鐘差建模及其在精密單點定位中的應用

張小紅，陳興漢，郭 斐

武漢大學測繪學院，湖北 武漢430079

1 引 言

近年來，精密單點定位（precise point positioning，PPP）技術的快速發展使其成為獲取高精度測站坐標的一種重要技術手段，已被廣泛應用于地球動力學、電離層和對流層延遲估計等諸多地學研究。但是，受衛星定位幾何條件的限制，GNSS精密定位獲得的測站高程精度明顯次于平面位置精度（1.5～2倍）。因此，在利用GNSS跟蹤站的觀測資料計算地殼垂直運動量時，一些微弱的地學信號往往被淹沒在測站坐標時間系列的噪聲中而無法準確提?。?－2］。造成 GNSS定位高程精度偏低的主要原因是地球遮擋影響，接收機捕獲的衛星信號僅來自地平線以上的可視衛星，這種幾何上非對稱性的觀測條件使得接收機鐘差、天頂對流層延遲以及測站高程參數之間存在顯著的數學相關性（瞬時相關性高達80%～90%），這種強相關性使得高程方向的定位精度偏低［3－4］。盡管引入低高度角衛星的觀測數據能夠在一定程度上削弱接收機鐘差、測站高程以及天頂對流層延遲三者之間的相關性，但是由于低高度角衛星受到的系統誤差（大氣殘余誤差、多路徑效應）和觀測噪聲（低信噪比）明顯較大，反而可能降低水平分量的定位精度。解決這一問題的一種有效途徑是引入外部約束信息。假設能夠使用一些高頻穩度的接收機鐘，如原子鐘或超穩石英鐘，充分利用其短期頻穩的約束條件，對接收機鐘差進行建模與預報，則有望顯著減少接收機鐘差參數個數，降低上述3種參數的相關性，從而更加準確地分離（估計）出測站高程與天頂對流層延遲。

長期以來，國內外學者對時鐘性能評估與鐘差建模方面的研究主要是針對衛星端高性能的原子鐘，并取得了豐碩的研究成果［5－10］。針對地面或星載GNSS接收機端的鐘差建模及其應用研究相對較少，尚處于起步階段。近年來，隨著IGS跟蹤站網硬件設備的不斷升級，目前已有超過130個跟蹤站（接收機端）配置了高準確度和高頻穩度的氫原子鐘、銣原子鐘或銫原子鐘［11］；此外，一些科學試驗衛星如GRACE重力衛星上也搭載了超穩振蕩器［12］，其短期（1～1000s）頻穩度高達1×10－13～3×10－13（文獻［13］）。這些高穩定度的接收機鐘使得地面或星載接收機鐘差建模成為可能。

因此，本文在分析評價當前IGS跟蹤站的幾類原子鐘性能的基礎上，擬重點研究并分析鐘差建模方法在精密單點定位中的應用。本文第2節將利用Allan方差法簡要評估現有IGS幾類原子鐘的穩定性能，第3節給出了本文鐘差建模的二次多項式模型，第4節重點討論了鐘差建模在精密度單點定位中的應用及效果。

2 接收機原子鐘的性能評估

鐘差性能（頻穩度）評估是接收機鐘差建模與預報的前提，目前常用的GNSS接收機鐘類型有石英鐘、銣原子鐘、銫原子鐘、氫原子鐘等［14］，但并非所有類型的接收機鐘都適合建模，特別是對于一些頻穩度較差的石英鐘，即便采用復雜的鐘差模型也難以準確描述時鐘的運行特性。此外，溫度變化、空間環境（地面和星載）差異等因素也會對接收機鐘的穩定性造成影響。表1給出了當前IGS跟蹤站原子鐘配置統計［15］，圖1為3類原子鐘的全球分布圖。

表1 全球IGS跟蹤站配備高穩原子鐘數目統計情況Tab.1 Number of atomic clocks in IGS tracking stations

圖1 配備高穩原子鐘的IGS跟蹤站全球分布圖Fig.1 IGS tracking stations with high stability atomic clocks

Allan方差（或標準差）是目前最常用的時域頻率穩定性分析方法［16］，它不僅可以用于計算接收機鐘的穩定度，還可用于識別接收機鐘的噪聲類型、計算噪聲水平系數。為了克服傳統的Allan方差無法識別調相白噪聲和調相閃爍噪聲這一缺陷［17］，本文采用修正Allan方差來表征接收機鐘的時域穩定度［18］，其計算公式為

式中，N為鐘差采樣數；τ為鐘差取樣（平滑）間隔；xi為歷元i對應的鐘差值；n為鐘差平滑因子，一般取為。

噪聲類型識別方面，通過修正Allan方差雙對數圖的斜率來區分噪聲過程，相應給出了振蕩器噪聲的修正Allan方差在雙對數圖上的表現形式，如圖2所示。從圖中可以得到，調相白噪聲（WPM）和調相閃爍噪聲（FPM）的斜率是－1，調頻白噪聲（WFM）的斜率是－0.5，調頻閃爍噪聲（FFM）的斜率是0，調頻隨機游走噪聲（RWFM）的斜率是0.5。

圖2 振蕩器噪聲類型的修正Allan方差Fig.2 The modified Allan variance of oscillator noise

為了分析接收機鐘的穩定性，首先采用精密單點定位獲得了部分IGS跟蹤站30s間隔的接收機鐘差，然后對接收機鐘差序列采用修正Allan方差計算得到部分IGS跟蹤站接收機鐘（H：氫原子鐘；Cs：銫原子鐘；Rb：銣原子鐘；QUARTZ：石英鐘）的半日穩定度。如圖3所示，其中橫軸代表時間間隔；縱軸代表對應的Allan標準方差數值（數值越小，穩定度越高）；兩條黑色虛線（GPS code clock，GPS phase clock）分別代表無電離層組合偽距和載波相位觀測值的等效噪聲下界［19］。

從圖3中可以看出，當采樣間隔在30～2000s左右時，4款接收機鐘的噪聲類型主要表現為調頻白噪聲，當采樣間隔取至2000s甚至更長時，調頻閃爍噪聲占主導地位。通過比較4種類型的接收機鐘差在雙對數圖的縱坐標值，可以明顯地看出氫原子鐘的穩定度比銫原子鐘要好，銫原子鐘比銣原子鐘的穩定度高，石英鐘的穩定度最差。分析圖3中無電離層組合偽距和載波相位觀測值的噪聲下界不難發現，石英鐘的頻穩度較差，即使在非常短的時間內也無法滿足載波相位水平的鐘差建模要求（即模型預報誤差大于觀測噪聲）；銫原子鐘具有較好的短期穩定度，但是當平滑間隔在超過300s時，其預報誤差也將超出載波相位觀測值的噪聲水平；而氫原子鐘的穩定度最高，可以滿足1～2h以內的鐘差建模精度要求。

圖3 4款接收機鐘差的修正Allan方差Fig.3 The modified Allan variance of four types of receiver clock offset

根據Allan方差原理，鐘差隨機項中不同噪聲分量引起的模型預報誤差RMSx（τ）存在所示的對應關系，見表2［16］。給定預報誤差的閾值RMSmax，利用表2中的模型預報誤差計算公式（RMSx（τ）≤RMSmax），即可以確定鐘差建模（模型參數）的有效時長。以圖3中氫原子鐘差建模為例，設定預報誤差的閾值RMSmax為0.006m，平滑間隔2000s以內的噪聲類型表現為調頻白噪聲，其預報誤差的計算公式為

式中，RMS（τ）為預報誤差；τ為預報的時間間隔；σy（τ）為平滑間隔為τ時的修正Allan標準差；c為光速。

表2 不同噪聲類型引起的模型預報誤差Tab.2 The prediction error of different types of noise

當τ取最大值2000s時，由圖3知σy（τ）為10－14，根據式（2）計算得到的預報誤差為0.006m，滿足RMS（τ）≤RMSmax的條件，可以在30～2000 s內任取采樣間隔進行鐘差建模。當采樣間隔超過2000s時，調頻閃爍噪聲占主導地位，調頻閃爍噪聲的預報誤差計算公式為

由圖3知，采樣間隔為2000～8000s以內的σy（τ）基本在10－14左右，根據式（3）算得的τ為1665s。綜合以上分析，得到跟蹤站mac2鐘差建模的有效時長取為2000s。

3 接收機鐘差模型

接收機鐘差可以用確定性變化分量和隨機性變化分量來描述［8］，即

式中，右邊前3項為鐘的確定性時間分量；a0、a1和a2依次代表接收機的鐘差、鐘速與鐘漂；εx（t）為接收機鐘差的隨機變化分量。接收機鐘差的系統性變化部分（確定性分量）可采用線性方程或二階多項式等確定性函數模型來表達［20］，而其隨機性變化分量只能從統計意義上來分析，其統計特性由噪聲冪律譜模型（式5）確定［21］。

瞬時時間偏差x（t）的功率譜密度可表示為

式中，fα（α＝－4，－3，－2，－1，0）代表5種噪聲類型的傅氏頻率；hβ（β＝－2，－1，0，1，2）為噪聲強度系數［22］。對于給定的原子鐘，隨著取樣時間的變化將表現出不同的隨機噪聲分量，但在較短的時間內，一般只有1～2種噪聲起主導作用。

4 鐘差建模在精密單點中的應用

4.1 附有鐘差約束的精密單點定位方法

在傳統的精密單點定位［23］中，采用Kalman濾波作為參數估計器，一般將接收機鐘差視作獨立的白噪聲（過程噪聲非常大），忽視了鐘差參數之間可能存在的短期相關性。因此，本文采用鐘差－鐘速二維狀態模型描述接收機鐘的動態過程［24］，其狀態一步預測方程可表示為

式中，xp和xf分別代表接收機鐘差、鐘速；Δt為歷元間隔；wp和wf為狀態的過程噪聲。

根據誤差傳播定律［25］，得到狀態一步預測的協因數矩陣Qx，k為

式中，Φ為狀態轉移矩陣；Qw，k為系統過程噪聲，取決于鐘的穩定度，可根據Allan方差或譜密度系數確定［26］

式中，Sp、Sf分別代表引起時差和頻差的噪聲譜振幅。

采用Kalman濾波算法即可獲得所有歷元的遞推解，為了確保所有歷元濾波解的精度和可靠性，本文采用雙向平滑濾波算法［27］對其結果進行優化。設前向濾波和后向濾波的解分別為Xk，f、Xk，b，對應的估值協因數矩陣為Qk，f、Qk，b，則平滑濾波的估值及其協因數矩陣為

4.2 試驗結果與分析

利用2010年4月18日WTRZ跟蹤站（配備氫原子鐘）的觀測數據（30s采樣率）和CODE分析中心提供的精密星歷和精密鐘差產品，依次采用以下兩種方案進行動態PPP解算：

方案1：采用傳統的逐歷元估計一維鐘差參數的方法，過程噪聲設置為3×105；

方案2：采用本文的鐘差建模方法估計鐘差－鐘速二維狀態［28］，過程噪聲由譜密度系數（本文根據氫原子鐘的特性設置了h0＝1×10－24；h－1＝4×10－29）確定。

基于上述兩種方案獲得的測站高程分量偏差、接收機鐘差參數以及測站高程分量與接收機鐘差參數之間的相關系數時序，如圖4—6所示。利用PPP逐歷元解算的方差協方差陣中的相關數據，計算得到各個歷元的相關系數ρ

式中，cov（xh，xt）代表高程分量與接收機鐘差參數的協方差；σxh代表高程分量的方差；σxt代表接收機鐘差參數的方差。

類似的，根據2010年4月18日WAB2跟蹤站（配備氫原子鐘）的觀測數據（30s采樣率）和CODE分析中心提供的精密星歷和精密鐘差產品采用上述兩種方案進行動態PPP解算，得到測站的高程分量偏差、高程分量與接收機鐘差參數之間相關系數如圖7和圖8所示。

圖4 兩種方案解算的高程方向偏差（WTRZ）Fig.4 Time series of height errors（WTRZ）

圖5 兩種方案解算的接收機鐘差參數（WTRZ）Fig.5 Time series of receiver clock errors（WTRZ）

圖6 高程分量與接收機鐘差參數的相關系數時序（WTRZ）Fig.6 Correlation coefficients between height component and receiver clock offset（WTRZ）

圖7 兩種方案解算的高程方向偏差（WAB2）Fig.7 Time series of height errors（WAB2）

圖8 高程分量與接收機鐘差參數的相關系數時序（WAB2）Fig.8 Correlation coefficient between height component and receiver clock offset（WAB2）

本文采用正反向平滑濾波算法（后處理），而不是單向濾波處理，即后處理PPP正反向平滑濾波算法消除了PPP單向濾波的收斂過程。分析上述結果不難發現，采用將接收機鐘差視為白噪聲的鐘差逐歷元估計方案，由其得到的測站高程分量、天頂對流層延遲與接收機鐘差參數之間存在顯著的相關性（特別是測站高程分量與接收機鐘差參數之間的相關性高達80%～90%），由此導致動態PPP高程方向的定位精度偏低，且波動較大，這就使得一些微弱的地學信號往往被淹沒在測站坐標時間系列的噪聲中而無法準確拾取。方案二則充分利用了接收機鐘差參數之間的短期穩定性，通過鐘差建模在一定程度上削弱了測站高程分量、天頂對流層延遲與接收機鐘差參數之間數學相關性，進而改善了動態PPP高程方向的定位精度，相應的RMS提高了50%左右。

此外，本文還對比分析了上述兩種方案估計的天頂對流層延遲精度。利用2010年4月18日AMC2、TWTF、WAB 2、WDC3、WTZR 5個測站（配備氫原子鐘）的GPS觀測數據進行動態PPP解算獲得了各觀測站的天頂對流層延遲參數（ZPD），并與IGS分析中心提供的ZPD參考值進行比較，統計兩種方案解算的天頂對流層延遲參數的外符合精度，見表3。結果表明，基于接收機鐘差建模的精密單點定位對各測站的ZPD參數估計均有不同程度的改善，尤其是WTZR和TWTF站的ZPD精度提高了近20%。

表3 兩種方案解算的天頂對流層延遲參數的精度比較Tab.3 Comparisons of the accuracy of tropospheric delay

對比前文高程分量的改善程度，接收機鐘差建模的方法對ZPD參數估值精度的改善幅度不及測站高程分量。這主要是因為ZPD參數與接收機鐘差參數之間的相關性次于測站高程分量與接收機鐘差參數之間的相關性。圖9給出了接收機鐘差參數與測站高程分量、ZPD參數的相關系數，虛線代表接收機鐘差參數與測站高程分量之間的相關系數，實線代表接收機鐘差參數與ZPD參數之間的相關系數，4條黑色虛線是判斷相關性強弱程度的臨界，接收機鐘差參數與ZPD參數的相關系數的絕對值基本維持在0.3～0.5之間；接收機鐘差參數與測站高程分量之間的相關系數的絕對值卻大于0.5。

5 結 論

針對當前許多IGS跟蹤站均配置有高性能原子鐘的這一現狀，本文首先采用修正Allan方差評估了不同類型接收機鐘的短期頻穩度及鐘差建模的可行性，結果表明，不同類型接收機鐘的穩定性有所差異，在給定模型預報誤差限值條件下（預報誤差小于觀測噪聲水平），銫原子鐘僅能滿足較短時間（數分鐘）內的建模精度要求，氫原子鐘具有良好的短期穩定度，能夠滿足1～2h內的鐘差建模精度要求。實例說明本文方法能有效改善動態PPP高程方向的定位精度，對天頂對流層延遲參數的估計也有一定改善。

［1］CRYER J D，CHAN K S.Time Series Analysis：with Applications in R［M］.2nd ed.Berlin：Springer，2008.

［2］GU Guohua.Recent Progress in Researches on Crustal Movement through GNSS（GPS）Observations［J］.Recent Developments in World Seismology，2007，343（7）：9－15.（顧國華.GNSS（GPS）觀測研究地殼運動的新進展［J］.國際地震動態，2007（7）：9－15.）

［3］ROTHACHER M，BEUTLER G.The Role of GPS in the Study of Global Change［J］.Physics and Chemistry of the Earth，1998，23（9）：1029－1040.

［4］DACH R，BEUTLER G，HUGENTOBLER U.Time Transfer Using GPS Carrier Phase：Error Propagation and Results［J］.Journal of Geodesy，2003，77（1）：1－14.

［5］HESSELBARTH A，WANNINGER L.Short－term Stability of GNSS Satellite Clocks and Its Effects on Precise Point Positioning［C］∥Proceedings of ION GNSS 2008.Savannah，GA：［s.n.］：1885－1863.

［6］HUANG Guanwen，ZHANG Qin，WANG Jigang.Research on Estimation and Prediction of GPS Satellite Clock Error［J］.Journal of Geodesy and Geodynamics，2009，29（6）：118－122.（黃觀文，張勤，王繼剛.GPS衛星鐘差的估計與預報研究［J］.大地測量與地球動力學，2009，29（6）：118－122.）

［7］GUO Hairong，YANG Yuaxi，HE Haibo，et al.Determination of Covariance Matrix of Kalman Filter Used for Time Prediction of Atomic Clocks of Navigation Satellites［J］.Acta Geodaetica et Cartographica Sinica，2010，39（2）：147－149.（郭海榮，楊元喜，何海波，等.導航衛星原子鐘Kalman濾波中噪聲方差－協方差的確定［J］.測繪學報，2010，39（2）：147－149.）

［8］DELPORTE J，BOULANGER C，MERCIER F.Short－term Stability of GNSS On－board Clocks Using the Polynomial Method［C］∥European Frequency and Time Forum（EFTF）.［S.l.］：IEEE，2012，117－121.

［9］GONG Hang，YANG Wenke，LIU Zengjun，et al.Estimation Method of BDS On－board Clock Short－term Stability Combining Satellite Two－way with One－way Carrier Ranging［J］.Journal of National University of Defense Technology，2013，35（3）：158－163.（龔航，楊文可，劉增軍，等.衛星雙向與單向載波聯合的北斗星載鐘短穩評估方法［J］.國防科技大學學報，2013，35（3）：158－163.）

［10］HAUSCHILD A，MONTENBRUCK O，STEIGENBERGER P.Short－term Analysis of GNSS Clocks［J］.GPS Solutions，2013，17（3）：295－307.

［11］WANG K，ROTHACHER M.Stochastic Modeling of Highstability Ground Clocks in GPS Analysis［J］.Journal of Geodesy，2013，87（5）：427－437.

［12］DUNN C，BERTIGER W，BAR－SEVER Y，et al.Application Challenge－instrument of Grace－GPS Augments Gravity Measurements：Twin Satellites Trail Each Other in Earth Orbit.As They Pass over Contours in the Gravity Field，They First［J］.GPS World，2003，14（2）：16－29.

［13］WEINBACH U，SCHON S.GNSS Receiver Clock Modeling when Using High－precision Oscillators and Its Impact on PPP［J］.Advances in Space Research，2011，47（2）：229－238.

［14］ALLAN D W.Statistics of Atomic Frequency Standards［J］.Proceedings of the IEEE，1966，54（2）：221－230.

［15］CLKLOG.A Summary File of the Deployment History for GPS Receiver，Antenna，Frequency Standards，and Other Equipment at IGS Stations［EB/OL］.2010－01－20［2014－05－30］ftp：∥igsws.unavco.org/igscb/station/general/loghist.txt.

［16］ALLAN D W.Time and Frequency（Time－domain）Characterization，Estimation，and Prediction of Precision Clocks and Oscillators［J］.IEEE Transactions on Ultrasonics，Ferroelectrics and Frequency Control，1987，34（6）：647－654.

［17］ALLAN D W，WEISS M A.A Frequency－domain View of Time－domain Characterization of Clocks and Time and Frequency Distribution Systems［C］∥Proceedings of the 45th Annual Symposium on Frequency Control：667－678.

［18］ALLAN D W，HOWE D A，WALLS F L.Characterization of Clocks and Oscillators［M］.Gaithersburg：US Department of Commerce，National Institute of Standards and Technology，1990.

［19］WEINBACH U，SCH?N S.Improved GRACE Kinematic Orbit Determination Using GPS Receiver Clock Modeling［J］.GPS Solutions，2013，17（4）：511－520.

［20］GUO Hairong.Research on Theory and Method about Time－frequency Characteristic of Atomic Clocks of Navigation Satellites［D］.Zhengzhou：Information Engineering University，2006：30－33.（郭海榮.導航衛星原子鐘時頻特性分析理論與方法研究［D］.鄭州：解放軍信息工程大學，2006：30－33.）

［21］BROWN R G，HWANG P Y.Introduction to Random Signals and Applied Kalman Filtering：with MATLAB Exercise and Solutions［J］.Wiley Global Education，2012（1）：1－10.

［22］BARNES J A，CHI A R，CUTLER L S，et al.Characterization of Frequency Stability［J］.IEEE Transactions on Instrumentation and Measurement，1971，1001（2）：105－120.

［23］KOUBA J.A Guide to Using International GNSS Service（IGS）Products［J］.International GNSS，2009，13（8）：1－10.

［24］STRANG G，BORRE K.Linear Algebra，Geodesy，and GPS［M］.Siam：Wellesley－Cambridge Press，1997.

［25］ZUMBERGE J F，HEFLIN M B，JEFFERSON D C，et al.Precise Point Positioning for the Efficient and Robust Analysis of GPS Data from Large Networks［J］.Journal of Geophysical Research，1997，102（B3）：5005－5017.

［26］HERRING T A，DAVIS J L，SHAPIRO I I.Geodesy by Radio Interferometry：The Application of Kalman Filtering to the Analysis of Very Long Baseline Interferometry Data［J］.Journal of Geophysical Research，1990，95（B8）：12561－12581.

［27］GELB A.Applied Optimal Estimation［M］.Cambridge：MIT Press，1974.

［28］SEEBER G.Satellite Geodesy：Foundations，Methods and Applications［M］.2nd ed.New York：Walter de Gruyter，2003.