一種改進的高精度頻率源鐘差仿真方法

徐 超，李 博，劉軍良，吳華兵，胡永輝

(1.中國科學院國家授時中心，西安 710600；2.中國科學院時間頻率基準重點實驗室，西安 710600；3.中國科學院大學，北京 100190)

一種改進的高精度頻率源鐘差仿真方法

徐 超1,2,3，李 博1,2,3，劉軍良1,2，吳華兵1,2，胡永輝1,2

(1.中國科學院國家授時中心，西安 710600；2.中國科學院時間頻率基準重點實驗室，西安 710600；3.中國科學院大學，北京 100190)

提出了一種適用于高精度頻率源仿真的鐘差仿真方法。從頻率源的噪聲特性及噪聲冪率譜模型入手，利用[-1,1]范圍內(nèi)服從均勻分布的隨機變量模擬各類噪聲的特性，并給出了各類噪聲模擬產(chǎn)生的數(shù)學表達式。在此基礎上，提出了基于牛頓迭代法求解噪聲參數(shù)的新方法，避免了噪聲參數(shù)擬合過程中進行線性近似帶來的誤差，并進行了試驗驗證。試驗結果表明，所提出的鐘差模擬方法適用于不同類型的原子鐘以及高穩(wěn)晶振的鐘差模擬仿真，模擬產(chǎn)生的鐘差數(shù)據(jù)的阿倫偏差與實測值非常接近，能很好地反映各種噪聲對鐘差數(shù)據(jù)的影響，可以滿足高精度頻率源模擬仿真的需求。

原子鐘；高穩(wěn)晶振；鐘差；冪率譜噪聲；牛頓迭代

0 引 言

高精度頻率源在時頻實驗室、導航、通信、深空探測等許多領域中應用越來越廣泛[1-3]，其鐘差仿真能夠為系統(tǒng)設計方案的評估、算法的仿真和應用的開發(fā)等提供可靠的驗證手段，因此也越來越受到系統(tǒng)設計者的重視。

對振蕩器噪聲的建模及仿真是鐘差仿真的基礎，文獻[4]對閃爍噪聲的特性進行了深入的研究，提出了使用白噪聲模擬產(chǎn)生閃爍噪聲的方法；文獻[5-7]分別對幾種常用的噪聲仿真模型進行了討論，并給出了仿真結果；文獻[8]提出了一種仿真冪律譜噪聲的新方法，其基本思路是將冪律譜噪聲x(t)看成是白噪聲w(t)通過一個線性系統(tǒng)的輸出，并給出了在連續(xù)域和離散域求解系統(tǒng)沖擊響應h(t)的方法；文獻[9]提出了一種利用服從均勻分布的隨機變量模擬冪律譜噪聲的簡便方法，并給出了時域表達式，文獻[10]對該方法進行了改進，加入了歸一化系數(shù)，使產(chǎn)生的各種噪聲的阿倫方差保持不變，并給出了使用線性最小二乘法求解噪聲參數(shù)的方法。

噪聲參數(shù)的選取對于最后仿真結果有著較大影響，不準確的噪聲參數(shù)會使得模擬產(chǎn)生的鐘差數(shù)據(jù)的阿倫方差與實測的阿倫方差產(chǎn)生較大偏離。本文提出了一種基于牛頓迭代法計算噪聲系數(shù)的方法，能準確的模擬頻率源的噪聲特性，與文獻[10]的方法相比，該方法模擬產(chǎn)生的鐘差數(shù)據(jù)的阿倫方差與實測值具有更好的符合度。

1 鐘差仿真建模

1.1鐘差信號的表示

鐘差指的是由頻率源瞬時相位變化引起的相對于參考系理想時間的相對時間偏差[11]。對于一個頻率源，瞬時頻率相對于標稱頻率的偏差y(t)定義為

(1)

式中：v0為標稱頻率，v(t)為輸出信號的瞬時頻率，y(t)表示輸出信號的瞬時頻率相對于標稱頻率的歸一化瞬時頻率偏差。相對頻率偏差y(t)在時域上的積分定義為鐘差x(t)：

(2)

(3)

1.2頻率源噪聲模型

對于一個頻率源而言，鐘差信號的不穩(wěn)定主要源于頻率源的物理特性和外部環(huán)境所引起的各種隨機噪聲。冪率譜噪聲模型是使用最廣泛的頻率源噪聲模型，它的有效性已經(jīng)被許多實驗結果所證實，并為國際學術界所普遍接受[11]。

冪率譜噪聲模型是用5種相互獨立的隨機噪聲來描述振蕩器的頻率波動過程，并通過輸出頻率的功率譜密度函數(shù)在對數(shù)域的斜率來區(qū)分各種不同的噪聲類型，功率譜密度函數(shù)表示為[12-13]

(4)

式中：hα為常數(shù)，α=-2,-1,0,1,2分別對應隨機游走噪聲、頻率閃爍噪聲、頻率白噪聲、相位閃爍噪聲、相位白噪聲5種不同的噪聲類型。

對于高精度頻率源鐘差模擬仿真而言，往往更關注這些噪聲在時間x(t)上的表現(xiàn)。由于時間是頻率的積分，x(t)的功率譜密度函數(shù)表示為：

(5)

將式(4)代入式(5)，可得

(6)

式中:功率譜密度函數(shù)Sx(f)在對數(shù)域的表示如圖1所示。

圖1 不同噪聲的功率譜密度Fig.1 Power spectral density of typical clock noise

如圖1所示，不同的曲線斜率可以區(qū)分不同類型的隨機噪聲。隨機游走噪聲是最低頻的噪聲，其對長期頻率穩(wěn)定度的影響較大，而從短期來看，相位白噪聲則對頻率穩(wěn)定度的影響最大。

阿倫方差是常用的頻率穩(wěn)定度的時域表征，定義為[13]

(7)

式中：τ表示采樣時間，N表示以τ為采樣間隔的樣本數(shù)。在實際應用中，一般采用阿倫方差的平方根來表征，即阿倫偏差。

阿倫方差與譜密度Sy(f)之間存在轉(zhuǎn)換關系，其表達式為

(8)

在實際測量過程中，系統(tǒng)輸出一般需要經(jīng)過一個低通濾波器，進入低通濾波器的頻率超過濾波器的截止頻率fh時便被抑制，此時表達式為

(9)

其中,低通濾波器的截止頻率fh滿足：

(10)

結合式(4)和式(9)可以得到冪率譜噪聲的阿倫方差，結果如表1所示，對于鐘差信號的仿真也是基于表1內(nèi)容展開的。

表1 各種噪聲功率譜密度與阿倫方差換算表Table 1 Conversion table of noise power spectral density and Allan variance

1.3頻率源噪聲模擬方法

一般情況，在進行鐘差數(shù)據(jù)模擬產(chǎn)生時，可以忽略相位閃爍噪聲[9-10]。本小節(jié)主要對其他四種噪聲進行分析。

令

作為四種噪聲的噪聲參數(shù)，將上述噪聲參數(shù)代入表1，各種噪聲與阿倫方差之間的函數(shù)關系如表2所示。

使用服從均勻分布的隨機變量能很好的模擬各種噪聲的特性，各種噪聲引起的相對頻率偏差的表達式為[9]

表2 各種噪聲與阿倫方差Table 2 Allan variance of types of clock noise

(11)

(12)

(13)

式中：D(i)為[-1,1]范圍內(nèi)的服從均勻分布的隨機變量，i=1,2,…,N，N為仿真數(shù)據(jù)點數(shù)，τs為采樣周期。

使用式(11)～(13)可以很容易的模擬產(chǎn)生相位白噪聲、頻率白噪聲和隨機游走噪聲的相對頻率偏差。然而，頻率閃爍噪聲的模擬產(chǎn)生則比較不容易，文獻[14]中提出了模擬產(chǎn)生頻率閃爍噪聲的方法，取得了很好的效果。其表達式為

(14)

在使用上述方法產(chǎn)生各種噪聲引起的相對頻率偏差數(shù)據(jù)時，必須加入一個歸一化常量k，如表3[10]所示，以使得產(chǎn)生的各種噪聲數(shù)據(jù)的阿倫方差保持不變。

表3 歸一化系數(shù)Table 3 Normalized coefficient

按照上述方法，能夠模擬產(chǎn)生各類噪聲的相對頻率偏差，再將式(11)～(14)代入式(2)，便可得到鐘差的模擬數(shù)據(jù)，鐘差數(shù)據(jù)產(chǎn)生的表達式為：

x(i)=x(i-1)+[ywp(i)+

ywf(i)+yff(i)+yrw(i)]τs

(15)

2 基于牛頓迭代法的噪聲參數(shù)擬合方法

文獻[10]提出了一種確定噪聲參數(shù)的方法。利用對頻率源的實際測量得到的鐘差數(shù)據(jù)，計算出多個阿倫偏差，再利用線性最小二乘方法及非負的約束條件，可以計算得到噪聲參數(shù)，其表達式為：

(16)

式中：fh=1 Hz，σy為阿倫偏差。

采用上述方法進行噪聲參數(shù)確定時，式(16)的各個等式是對阿倫偏差表達式(17)的線性近似，其阿倫偏差與理論值存在一定的偏差，如圖2所示，這一偏差將對噪聲參數(shù)的擬合造成較大的影響。

(17)

圖2 某典型氫原子鐘的阿倫偏差Fig.2 Allan standard deviation of a typical hydrogen atomic clock

在實際應用過程中，需要使得模擬產(chǎn)生的鐘差數(shù)據(jù)符合某特定頻率源的統(tǒng)計特性。因此，將該頻率源的阿倫偏差作為輸入，將式(16)按照式(17) 進行修正，得到如下超定非線性方程組。

(18)

牛頓迭代法是求解非線性方程組的常用有效方法之一，可以用于求解上述超定非線性方程組的最優(yōu)解A=[Awp,Awf,Aff,Arw][15]。

(19)

超定非線性方程組的一般形式如式(19)所示，其中,fi(x1,x2,…,xn)是定義在區(qū)域D?Rn上的實值函數(shù)，m>n。按照牛頓迭代法求得方程組迭代k+1次后的解Xk+1為

(20)

(21)

在式(20)中，G(Xk)被加上了下標，是因為每一次更新Xk+1以后，F(xiàn)(X)的Jacobi矩陣都要重新計算。更新結束的條件是通過判斷X是否收斂，即當修正量ΔXk滿足

ΔXk<ΔXth

(22)

其中,ΔXth是預先設定的一個閾值，當ΔXk小于該閾值時，就認為可以停止更新了。X的收斂只是意味著我們已經(jīng)找到了使最小二乘的代價函數(shù)最小的解，并不是說代價函數(shù)已經(jīng)趨于0。

在使用牛頓迭代法解非線性方程組時，如果迭代矩陣G(Xk)出現(xiàn)奇異或病態(tài)情況，迭代式(20)往往無法進行或者很難奏效。因此，為了解決迭代過程中G(Xk)出現(xiàn)奇異或病態(tài)的問題，可以采用阻尼牛頓法對牛頓迭代法進行修正[15]。

阻尼牛頓法是在牛頓迭代法的基礎上增加一個阻尼因子α∈(0,1]，使得

Xk+1=Xk+αp

(23)

(24)

圖3 阻尼牛頓法求解非線性方程組流程圖Fig.3 The flow chart of solve nonliner equation using damped Newton iteration

3 試驗結果

3.1氫原子鐘試驗結果

試驗利用國家授時中心守時實驗室的兩臺氫原子鐘MHM2010，采用互比法[18]測量其頻率穩(wěn)定度。利用本文中的方法模擬產(chǎn)生鐘差數(shù)據(jù)，計算其阿倫偏差，以MHM2010的實測阿倫偏差數(shù)據(jù)為參考，比較兩種鐘差仿真方法產(chǎn)生的鐘差數(shù)據(jù)的阿倫偏差與參考值的偏差，試驗結果見圖4和表4。

從圖4和表4可以看出，與文獻[10]的采用最小二乘法確定噪聲參數(shù)的鐘差仿真方法相比，采用本文提出的基于牛頓迭代法的鐘差仿真方法模擬產(chǎn)生的鐘差數(shù)據(jù)的阿倫偏差與實測值之間的偏差更小，能更好地反映各種噪聲對鐘差數(shù)據(jù)的影響。

圖4 某臺氫鐘MHM2010模擬鐘差數(shù)據(jù)阿倫偏差與實測值Fig.4 Allan deviation of a MHM2010 hydrogen atomic clock

τ/s實測值最小二乘法牛頓迭代法12×10-131.64×10-132.00×10-13103×10-142.03×10-142.94×10-141006×10-154.44×10-157.65×10-1510003×10-151.81×10-153.00×10-15100002.5×10-151.38×10-152.17×10-15864002×10-151.50×10-152.65×10-15

3.2銫原子鐘試驗結果

試驗采用以國家授時中心守時實驗室的主鐘信號為參考，利用高性能相噪與阿倫偏差分析儀5120A測得的某臺HP5071A銫原子鐘的頻率穩(wěn)定度實測數(shù)據(jù)。按照本文中的鐘差仿真方法模擬產(chǎn)生鐘差數(shù)據(jù)，計算其阿倫偏差，以HP5071A的實測阿倫偏差數(shù)據(jù)為參考，比較兩種鐘差仿真方法產(chǎn)生的鐘差數(shù)據(jù)的阿倫偏差與參考值的偏差，試驗結果見圖5和表5。

圖5 某臺銫鐘HP5071A模擬鐘差數(shù)據(jù)阿倫偏差與實測值Fig.5 Allan deviation of a HP5071A cesium atomic clock

τ/s實測值最小二乘法牛頓迭代法103.00×10-122.53×10-123.01×10-12202.03×10-121.77×10-122.03×10-12401.38×10-121.25×10-121.39×10-121008.7×10-137.88×10-138.67×10-132005.99×10-135.59×10-136.10×10-134004.06×10-133.94×10-134.31×10-1310002.7×10-132.51×10-132.73×10-1320001.93×10-131.80×10-131.94×10-1340001.12×10-131.32×10-131.39×10-13100005.2×10-141.0×10-139.06×10-14

與圖4氫原子鐘的試驗結果相比，兩種方法仿真產(chǎn)生的鐘差數(shù)據(jù)的阿倫偏差與實測值的偏差相差不大，本文提出的鐘差仿真方法的擬合精度只略優(yōu)于文獻[10]的采用最小二乘法確定噪聲參數(shù)的鐘差仿真方法。這是由于τ≤10000 s區(qū)間內(nèi)，銫原子鐘主要受頻率白噪聲的影響[19]，此時其他類型的噪聲對頻率穩(wěn)定度的影響較小，所以式(16)在對式(17)線性近似時引入的誤差較小。

3.3銣原子鐘試驗結果

與第3.2節(jié)銫原子鐘試驗方法類似，以同樣方法測得的某臺PRS10銣原子鐘的頻率穩(wěn)定度實測數(shù)據(jù)為參考，按照兩種鐘差仿真方法分別產(chǎn)生鐘差數(shù)據(jù)，計算并比較其阿倫偏差與參考值的偏差，結果見圖6和表6。

圖6 某臺銣鐘PRS10模擬鐘差數(shù)據(jù)阿倫偏差與實測值Fig.6 Allan deviation of a PRS10 rubidium atomic clock

圖6和表6的試驗結果表明，銣原子鐘的試驗結果與氫原子鐘的試驗結果類似。本文提出的鐘差仿真方法較之文獻[10]提出的采用最小二乘法確定噪聲參數(shù)的鐘差仿真方法有較為明顯的改善，模擬產(chǎn)生的鐘差數(shù)據(jù)的阿倫偏差與實測值更接近，能更好地體現(xiàn)各種噪聲對鐘差數(shù)據(jù)的影響。

表6 某臺銣鐘PRS10模擬鐘差數(shù)據(jù)阿倫偏差與實測值Table 6 Allan deviation of a PRS10 rubidium atomic clock

3.4高穩(wěn)晶振試驗結果

高穩(wěn)晶振OSA8607的頻率穩(wěn)定度數(shù)據(jù)來源與第3.2節(jié)、第3.3節(jié)相同，均為實際測量數(shù)據(jù)。采用與前文同樣的試驗方法，比較兩種方法產(chǎn)生的鐘差數(shù)據(jù)的阿倫偏差與實測值的偏差，結果見圖7和表7。

圖7 某高穩(wěn)晶振OSA8607模擬鐘差數(shù)據(jù)阿倫偏差與實測值Fig.7 Allan deviation of a OSA8607 high stability OCXO

τ/s實測值最小二乘法牛頓迭代法0.11.6×10-131.83×10-131.87×10-1317×10-145.87×10-147.6×10-14107.1×10-143.14×10-146.62×10-141008×10-147.07×10-149.98×10-1410002.6×10-132.19×10-132.60×10-13

采用最小二乘法模擬產(chǎn)生的鐘差數(shù)據(jù)的阿倫偏差與實測的頻率源的阿倫偏差在τ=1～100 s區(qū)間存在較大偏差，已經(jīng)不能反映高穩(wěn)晶振實際的頻率穩(wěn)定度特性，而本文提出的原子鐘鐘差模擬方法模擬產(chǎn)生的鐘差數(shù)據(jù)的阿倫偏差與實測值偏差較小，在統(tǒng)計學意義上能很好的反映各種噪聲對鐘差數(shù)據(jù)的影響。

綜上所述，與文獻[10]提出的采用最小二乘法確定噪聲參數(shù)的鐘差模擬方法相比，本文所提出的鐘差仿真方法在對噪聲參數(shù)的擬合建模過程中，避免了由于對式(17)進行線性近似而帶來的誤差，更加符合其頻率穩(wěn)定度的統(tǒng)計特性。模擬產(chǎn)生的鐘差數(shù)據(jù)的阿倫偏差與實測值的更加接近，能更好的反映各種噪聲對鐘差數(shù)據(jù)的影響，適用于不同類型的原子鐘和高穩(wěn)晶振的鐘差模擬仿真。

4 結 論

本文從頻率源的噪聲特性及噪聲冪率譜模型入手，提出了一種準確且簡便的模擬高精度頻率源鐘差信號的方法，模擬產(chǎn)生鐘差數(shù)據(jù)的阿倫偏差與給定的頻率源的阿倫偏差幾乎一致。該方法允許在不失實際隨機行為的情況下，對高精度頻率源的噪聲行為進行深度的分析，能夠方便的應用于時頻實驗室、衛(wèi)星導航系統(tǒng)、高精度通信系統(tǒng)、深空探測等許多領域的算法驗證和應用開發(fā)。

[1] 唐升,劉婭,李孝輝. 星載原子鐘自主完好性監(jiān)測方法研究[J]. 宇航學報, 2013, 34(1):39-45.[Tang Sheng,Liu Ya,Li Xiao-hui. A study on onboard satellite atomic clock autonomous integrity monitoring[J]. Journal of Astronautics, 2013, 34(1):39-45.]

[2] 張小紅,陳興漢,郭斐. 高性能原子鐘鐘差建模及其在精密單點定位中的應用[J]. 測繪學報, 2015, 44(4):392-398. [Zhang Xiao-hong, Chen Xing-han, Guo Fei. High performance atomic clock modeling and its application in precise point positioning[J]. Acta Geodaetica et Cartographica Sinica, 2015, 44(4):392-398.]

[3] Sch?fer W, Devicente J. ESA′s frequency and timing systems for deep space operations and radio science investigations[C]. The 20th European Frequency and Time Forum, Braunschweig, Germany, March 27-30, 2006.

[4] Barnes J A, Jarvis S. Efficient numerical and analog modeling of flicker noise processes[R]. Washington DC: U.S. Government Printing Office, June 1971.

[5] Kartaschoff P. Computer simulation of the conventional clock model[J]. IEEE Transactions on Instrumentation and Measurement, 1979, 28(3): 193-197.

[6] Radeka V. 1/|f| noise in physical measurements[J]. IEEE Transactions on Nuclear Science, 1969, 16(5): 17-35.

[7] Riggs T L, Phillips C L. Modeling continuous noise sources in digital simulations[J]. Simulation, 1987, 48(1): 11-18.

[8] Kasdin N J, Walter T. Discrete simulation of power law noise[C]. The 46th Frequency Control Symposium, Hershey, PA, USA, May 27-29, 1992.

[9] Harting A. Considering clock errors in numerical simulations[J]. IEEE Transactions on Instrumentation and Measurement, 1996, 45(3): 715-720.

[10] Diez J, D′Angelo P, Fernández A. Clock errors simulation and characterisation[C]. ION GNSS 19th International Technical Meeting of Satellite Division,Fort Worth,TX Sep 26-29,2006.

[11] 漆貫榮. 時間科學基礎[M]. 北京：高等教育出版社, 2006.

[12] 翟造成, 張為群, 蔡勇, 等. 原子鐘基本原理與時頻測量技術[M].上海: 上?？茖W技術文獻出版社,2009.

[13] Barnes J A, Chi A R, Cutler L S, et al. Characterization of frequency stability[J]. IEEE Transactions on Instrumentation and Measurement, 1971, 20(2): 105-120.

[14] Barnes J A, Allan D W. A statistical model of flicker noise[J]. Proceedings of the IEEE, 1966, 54(2): 176-178.

[15] 陳淑銘,喬田田. 一個求解非線性最小二乘問題的新方法[J]. 煙臺大學學報: 自然科學與工程版, 2004, 17(1): 14-22. [Chen Shu-ming, Qiao Tian-tian. A new method for solving problems of nonlinear least square [J]. Journal of Yantai University: Natural Science and Engineering Edition, 2004, 17(1): 14-22.]

[16] 張建軍,李春泉,張烈輝. 影響阻尼牛頓法收斂性的兩個重要參數(shù)[J]. 純粹數(shù)學與應用數(shù)學, 2012, 28(4): 433-439. [Zhang Jian-jun, Li Chun-quan, Zhang Lei-hui. The effect of two important parameter upon the convergence in damped newton method[J]. Pure and Applied Mathematics, 2012,28(4):433-439.]

[17] 朱小飛. 迭代法解非線性方程組[D]. 合肥：合肥工業(yè)大學, 2014. [Zhu Xiao-fei. Some iterative methods for solving system of nonlinear equations[D]. Hefei: Hefei University of Technology, 2014.]

[18] JJG 1004—2005, 氫原子頻率標準檢定規(guī)程[S].

[19] 張敏. 原子鐘噪聲類型和頻率穩(wěn)定度估計的自由度分析與探討[D]. 北京：中國科學院研究生院, 2008. [Zhang Ming. Analysis of the noise type of atomic clock and the degree of freedom of frequency stability[D]. Beijing: Graduate University of Chinese Academy of Sciences, 2008.]

AModifiedMethodforClockErrorSimulationofHighPrecisionFrequencySource

XU Chao1,2,3, LI Bo1,2,3, LIU Jun-liang1,2, WU Hua-bing1,2, HU Yong-hui1,2

(1.National Time Service Center, Chinese Academy of Sciences, Xi’an 710600, China;2.Key Laboratory of Time and Frequency Primary Standards, Chinese Academy of Sciences, Xi’an 710600, China;3.University of the Chinese Academy of Sciences, Beijing 100190, China)

An advanced approach to generating clock errors is presented that fits a high precision frequency source. Starting from the noise characteristics and noise power spectrum model, a new method for solving the noise parameters based on the Newton iterative method is proposed. Using these results as input, a sequence of clock errors are then obtained using the random variable subject to the uniform distributionin the range of [-1,1].Then, an experiment is carried out to test the performance of this simulation method. The Allan deviation of the simulation data is close to the reference value, and can reflect the effect of all kinds of noise on clock errors. It refers that the proposed simulation method is applicable for the clock error simulation of different types of atomic clocks and high stability crystal oscillators.

Atomic clock; High stability OCXO; Clock error; Power law noise; Newton iterative

TH714.1+4

A

1000-1328(2017)09- 0998- 08

10.3873/j.issn.1000-1328.2017.09.013

2017- 05- 11;

2017- 07- 17

中國科學院西部青年學者項目(XAB2015B13); 中國科學院西部青年學者資助項目(XAB2015B16)

徐超(1989-)，男，博士生，主要從事高精度時間頻率測量控制方法和技術研究。

通信地址：陜西省西安市臨潼區(qū)書院東路3號中國科學院國家授時中心(710600)

電話：(029)83890244

E-mail:xuchao1226@126.com