張小紅，陳興漢，郭 斐

武漢大學(xué)測繪學(xué)院，湖北 武漢430079

1 引 言

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

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

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

2 接收機(jī)原子鐘的性能評估

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

表1 全球IGS跟蹤站配備高穩(wěn)原子鐘數(shù)目統(tǒng)計情況Tab.1 Number of atomic clocks in IGS tracking stations

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

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

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

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

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

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

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

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

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

式中，RMS（τ）為預(yù)報誤差；τ為預(yù)報的時間間隔；σy（τ）為平滑間隔為τ時的修正Allan標(biāo)準(zhǔn)差；c為光速。

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

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

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

3 接收機(jī)鐘差模型

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

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

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

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

4 鐘差建模在精密單點中的應(yīng)用

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

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

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

根據(jù)誤差傳播定律［25］，得到狀態(tài)一步預(yù)測的協(xié)因數(shù)矩陣Qx，k為

式中，Φ為狀態(tài)轉(zhuǎn)移矩陣；Qw，k為系統(tǒng)過程噪聲，取決于鐘的穩(wěn)定度，可根據(jù)Allan方差或譜密度系數(shù)確定［26］

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

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

4.2 試驗結(jié)果與分析

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

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

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

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

式中，cov（xh，xt）代表高程分量與接收機(jī)鐘差參數(shù)的協(xié)方差；σxh代表高程分量的方差；σxt代表接收機(jī)鐘差參數(shù)的方差。

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

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

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

圖6 高程分量與接收機(jī)鐘差參數(shù)的相關(guān)系數(shù)時序（WTRZ）Fig.6 Correlation coefficients between height component and receiver clock offset（WTRZ）

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

圖8 高程分量與接收機(jī)鐘差參數(shù)的相關(guān)系數(shù)時序（WAB2）Fig.8 Correlation coefficient between height component and receiver clock offset（WAB2）

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

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

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

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

5 結(jié) 論

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

［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ìn)展［J］.國際地震動態(tài)，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衛(wèi)星鐘差的估計與預(yù)報研究［J］.大地測量與地球動力學(xué)，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.（郭海榮，楊元喜，何海波，等.導(dǎo)航衛(wèi)星原子鐘Kalman濾波中噪聲方差－協(xié)方差的確定［J］.測繪學(xué)報，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.（龔航，楊文可，劉增軍，等.衛(wèi)星雙向與單向載波聯(lián)合的北斗星載鐘短穩(wěn)評估方法［J］.國防科技大學(xué)學(xué)報，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，F(xiàn)requency 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，F(xiàn)erroelectrics 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ǎo)航衛(wèi)星原子鐘時頻特性分析理論與方法研究［D］.鄭州：解放軍信息工程大學(xué)，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.