吳 瓊滕云田黃大倫龍劍鋒

1）中國北京100081中國地震局地球物理研究所

2）中國北京100013中國計量科學(xué)研究院

絕對重力儀研制中一種新的自由落體軌跡重建算法＊

吳 瓊1），滕云田1）黃大倫2）龍劍鋒1）

1）中國北京100081中國地震局地球物理研究所

2）中國北京100013中國計量科學(xué)研究院

真空系統(tǒng)中落體做自由下落運動的軌跡重建算法，是進行擬合求解絕對重力測量結(jié)果的關(guān)鍵步驟之一.在自主研制小型高精度激光干涉絕對重力儀的過程中，基于目前高速發(fā)展的數(shù)字測量技術(shù)，提出一種自由落體軌跡重建的算法.其基本原理是：首先通過高速采樣，得到落體自由下落過程中產(chǎn)生的一系列數(shù)字化的激光干涉信號；其次通過數(shù)字處理算法提取干涉信號每次過零前后采樣點的時間坐標；然后按照相同的下落高度均分這些過零點，得到落體運動軌跡上的“時間－位移”坐標，從而重建落體軌跡.模擬實驗表明，該算法時間平均誤差在±0.01 ns范圍以內(nèi)，由此引入的計算測量重力加速度平均誤差為5.2×10－12m/s2，一次完整測量循環(huán)的時間為15 s，計算時間為2.5 s.該結(jié)果保證了高精度絕對重力儀研制中對測量速度和計算精度的要求.

絕對重力儀 自由落體 過零點計算 軌跡重建

Abstract：Rebuilding the trace of a free－fall body in the vacuum is one of the key steps in finding gravitational acceleration from the sampled fringe signals.This paper proposed an algorithm for rebuilding the trace during the independent process of developing small Laser Interference Absolute Gravimeter.The basic principles are：①Computing the time of the points before and after zero－crossing of a completely sampled fringe signal，which is generated by a freely falling body，and getting the time coordinates of zero－crossing points；②Dividing these time coordinates to get the Time－Displacement Coordinates in the bodytrace using the same displacement interval.In the experiment，we found that it needs 2.5 s to complete one calculation，the average timing accuracy is ±0.01 ns，the induced error is 5.2×10－12m/s2，and one complete measurement time is 15 s.The results satisfy the demand of the measurement speed and the accuracy in developing the high precision laser interference absolute gravimeter.

Key words：absolute gravimeter；free－fall body；calculation of zero－crossing；trace rebuilding

引言

絕對重力儀是一種高精度的重力加速度絕對值（10－8m/s2）測量設(shè)備，廣泛應(yīng)用于地震監(jiān)測（劉棟勛，劉文泰，1984；游澤霖，徐菊生，1984；Simon et al，2001），地?zé)豳Y源探測（Sugihara，2003），海平面測量（Vilibic，1997），相對重力儀的標定（Hinderer et al，1991），資源探測（張善法等，2009），以及潛艇的導(dǎo)航和導(dǎo)彈制導(dǎo)系統(tǒng)（程力等，2007；郭有光等，2003）等領(lǐng)域.一般絕對重力儀按照落體的運動方式可以分為上拋法和下落法（Sakuma，1984；Rothleitner et al，2009），g值的確定是通過追蹤計算位于真空室內(nèi)做自由下落（或上拋）運動的落體（或拋體）的運動軌跡得到的.算法的主要過程是首先獲得落體相對于參考點位置（參考棱鏡）信息的干涉帶信息，然后提取干涉帶上包含的落體運動的信息，通過最小二乘擬合獲得儀器歸算點處的重力加速度值（Durando et al，2002）.這個過程中落體運動信息的提取精度是保證結(jié)果精度的前提.2009年第8次國際絕對重力儀比對大會上，基于原子干涉的絕對重力儀第一次出現(xiàn)（de Angelis et al，2009）.這種儀器是通過磁光阱冷卻原子云團，從而提取原子團自由下落時感受的重力場信息.這種測量方式從原理上與激光干涉法測量的精度處于同一量級.

激光干涉絕對重力儀中，對干涉帶的處理可以分為硬件和軟件兩類.美國的JILA－g以及FG5系列重力儀均是基于硬件的過零觸發(fā)法提取干涉帶的過零點信息（Zumberge et al，1982；Van Camp et al，2003）.其主要過程是通過過零檢測器檢測干涉帶的過零點，然后將干涉帶過零時產(chǎn)生的每N個脈沖信號生成1個脈沖，送入與外部銣原子鐘同步的時間間隔分析儀中，從而測量每N個干涉帶所經(jīng)歷的時間間隔.這種方法主要是通過硬件的方式提取落體運動軌跡上的時間－位移對，然后利用最小二乘原理擬合落體運動軌跡，計算重力加速度.

隨著現(xiàn)代數(shù)字技術(shù)的高速發(fā)展，干涉帶的數(shù)字化處理與基于硬件的過零觸發(fā)處理電路相比，有著硬件系統(tǒng)和處理算法簡單、誤差小、模塊化的優(yōu)勢.意大利學(xué)者D’Agostino等提出了一種數(shù)字化干涉帶處理算法局部擬合算法（D’Agostino et al，2005，2008），通過數(shù)字處理算法提取干涉帶上落體感受到的重力作用信息.Bich等（2007）對該算法的測量誤差進行了詳細分析.該算法首先通過高速數(shù)字化儀對干涉帶信號進行采集，然后寫出一個帶有重力加速度變量的函數(shù)方程，通過總體最小二乘算法，與干涉帶信號直接擬合求解g值.為了減小誤差，將干涉信號分成若干段，每段內(nèi)利用函數(shù)方程擬合干涉波形，直接求解g值.由于每段包含了很短的采樣時間，故可以認為每段內(nèi)干涉信號的頻率恒定，這樣的近似的確可以減小調(diào)頻帶來的誤差.但實際落體下落過程中，干涉帶的頻率隨下落距離的加大而加快，局部擬合算法基于段內(nèi)干涉帶波形頻率固定的假設(shè)和對整個下落距離干涉帶擬合，無疑會增加算法的誤差和計算時間.

本文首先描述了筆者所在的研究團隊設(shè)計絕對重力儀的數(shù)據(jù)采集系統(tǒng)，然后介紹了該系統(tǒng)提出的一種全新的落體軌跡重建算法.該算法利用數(shù)字處理算法，提取數(shù)字化的干涉帶過零點信息，重建落體運動軌跡.通過模擬實驗，得到了這種算法的計算精度和所需時間.結(jié)果表明，這種算法在保證精度的前提下提高了運算速度，實現(xiàn)了對落體軌跡精確重建.最后給出了實驗樣機在2011年5月21日—6月2日的實驗結(jié)果.

1算法原理

在地球重力場作用下，真空系統(tǒng)中落體做自由下落運動時，反應(yīng)落體位置信息的干涉光束被光電接收器接收.光電接收器輸出電信號被高速數(shù)字化儀接收，轉(zhuǎn)變成數(shù)字信號后保存在數(shù)字化儀的內(nèi)存中.數(shù)字化儀的時間基準用外部銣原子鐘鎖定，如圖1所示.本系統(tǒng)采用的光源為基于非飽和蒸汽壓的碘穩(wěn)定633 nm He－Ne激光系統(tǒng)，工作時鎖定在I2分子飽和吸收譜線的f峰（錢進等，2008），頻率日長期穩(wěn)定性優(yōu)于5×10－11；光電接收器頻帶為4 kHz—80 MHz；銣原子鐘頻率重現(xiàn)性小于5×10－11，穩(wěn)定度小于1×10－11/s；高速數(shù)字化儀采樣頻率為62.5 MHz，時間基準通過鎖相環(huán)由外部銣原子鐘提供，板載內(nèi)存128 M，字長8位.

圖1 算法中的功能模塊連接示意圖Fig.1 Relationship between the function modules in the algorithm

落體相對于參考棱鏡位置信息的干涉帶波形可以認為是調(diào)頻的正弦波，其瞬時頻率隨著落體下落的速度的增加而不斷增加，下落的位移與干涉帶波形瞬時頻率的關(guān)系可以表示為

式中，λ為激光波長，本研究中的激光波長取633 nm；f為干涉帶波形的瞬時頻率；T為瞬時周期；s為落體自由下落的位移；g為重力加速度（1g＝9.8 m/s2）.

由此可以計算出，當(dāng)計算的落體下落距離為1—18cm時，干涉帶的頻率范圍為1.4—5.9 MHz.為了從得到的數(shù)字干涉帶波形中重建落體的下落軌跡，需要計算落體在某個時刻的下落距離.由于落體每下落半個激光波長，就會產(chǎn)生一條完整的干涉帶波形，同時數(shù)字化儀的時間基準來自于外部的銣原子鐘，因此可以通過計算一個完整的正弦波上的采樣點來得到落體運動軌跡上的時間－位移坐標，進而重建落體的下落軌跡.

因為干涉帶波形是一種調(diào)頻的正弦波形，同時又混入了各種背景噪聲，為了保證計算的精確度，可以通過計算信號斜率最大的過零點位置來重建落體的下落軌跡.在計算中，我們將落體的整個下落過程分成L段，每段中包含N個過零點，即每段對應(yīng)落體下落的距離為

式中，S為每段對應(yīng)落體自由下落的距離；N＝400，為每段中包含的干涉帶過零點次數(shù)；λ＝633 nm，為激光波長.由此可以計算出每段落體自由下落的位移約為0.06mm.

由于每段內(nèi)部我們僅需要統(tǒng)計干涉帶波形的過零點次數(shù)，所以不需要對每個過零點都精確計算其過零點時刻.因此，我們僅對每段的第一個過零點進行精確計算，具體計算過程如下：

1）首先設(shè)定一個噪聲上限閾值T，將絕對值大于T的采樣點組成一個新的時間序列Xnew，計算Xnew的過零點次數(shù)，然后提取每次過零點的前一個和后一個幅值大于T坐標，組成兩個新的時間序列Xind1和Xind2.圖2分別給出了過零點前（圓圈）和過零點后（星號）得到的Xind1和Xind2.Xind1的長度就是總過零點的個數(shù).

圖2 過零點算法結(jié)果（過零點前用圓圈表示，過零點后用星號表示）Fig.2 Result of the zero－crossing algorithm.Circles stand for the points before zero－crossing points，and stars denote the points after zero－crossing points

2）設(shè)定總段數(shù)L和每段內(nèi)過零點的個數(shù)N，由此確定每段落體的位移為S＝（N－1）/2 ×（λ/2），總位移為s＝S×L；

3）利用Xind1和Xind2計算每段第一個過零點的準確時刻.設(shè)干涉帶的采樣率為fs，對于第i段第一個過零點，該時刻的準確計算可以分為以下3種情況：

a 如果Xind2（i）與Xind1（i）之間無采樣點，則該時刻為

b 如果Xind2（i）與Xind1（i）之間僅有一個采樣點，在原始時間序列中的下標為k，則該時刻為

c 如果Xind2（i）與Xind1（i）之間多于一個采樣點，則對Xind2（i）與Xind1（i）之間的所有點進行置零，在Xind2（i）與Xind1（i）之間進行線性插值，設(shè)得到的插值曲線上最接近零的點下標為k，則該時刻為

2 模擬計算軌跡重建算法的精度

為了得到本文第一部分中描述的算法精度和效率，首先在雙核、4GB RAM的硬件環(huán)境測試.本算法采樣、計算總時間不超過15s，純計算時間不超過2.5s，如果考慮系統(tǒng)穩(wěn)定性要求而等待的時間，40s內(nèi)可以完成一次測量，每小時可以進行約90次測量.然后利用第一部分的計算過程1），通過不同頻率的正弦波來模擬落體在不同下落位置時產(chǎn)生的干涉帶信號，計算本算法在各個頻點處的誤差.

具體過程為：首先將銣原子鐘輸出的10MHz標準正弦波送入頻率綜合器，然后利用頻率綜合器輸出不同頻率的正弦波信號，通過集成有本算法的數(shù)字采集與處理系統(tǒng)分析計算這些正弦波信號，從而估計算法精度.這里用于實驗的頻率綜合器的輸出正弦波的時間誤差與外部銣原子鐘的時間不確定性一致.采集L段這樣的正弦波，每段內(nèi)有N個過零點，即（N－1）/2個完整的正弦波.設(shè)送入數(shù)字化儀的正弦波脈沖頻率為f，采樣頻率為fs，則每段正弦波對應(yīng)的標準時間可以表示為

通過第一部分的算法，可以計算每段正弦波起始點到終止點之間的采樣點個數(shù)為n，則這段波形對應(yīng)的時間可以表示為

則算法的時間誤差可以表示為

由此引入的絕對重力加速度值的計算誤差可以表示為

正弦波輸入數(shù)字化儀后，采集1000段，每段400個過零點，由式（8）計算得到1000段累積時間誤差、最大時間誤差、最小時間誤差和平均時間誤差.圖3顯示了落體從2cm即干涉帶瞬時頻率為1.98MHz開始，下落至18cm的過程中不同頻率點處算法的累積，最大和最小時間誤差，以及由此引入的測量誤差.

圖3 算法時間誤差（a）以及由此引入的測量誤差（b）估計Fig.3 Estimation of the algorithm time（a）and measurement（b）errors

通過利用式（1）和式（9）對算法時間誤差的估計發(fā)現(xiàn)，落體下落的高度在2cm之前，算法誤差較大；當(dāng)落體下落的高度超過2cm后，算法的累積時間誤差穩(wěn)定在±20ns之間；而當(dāng)落體下落高度大于6cm時，累積時間誤差穩(wěn)定在±10 ns之間.

由于實際干涉帶波形的頻率隨著下落高度的增加而增加，在相同的時間段內(nèi)，低頻段的過零點個數(shù)要遠遠小于高頻段，因此其過零點計算誤差對計算結(jié)果影響遠小于高頻段.因此在計算該算法在落體整個下落過程中誤差的綜合時，可以根據(jù)頻率設(shè)計一個加權(quán)系數(shù)，即

式中，fi是某個高度對應(yīng)的干涉帶信號瞬時頻率，wi為該頻率對應(yīng)的權(quán)系數(shù).表1給出了落體在1—18cm段、2—18cm段和6—18cm段內(nèi)算法對應(yīng)的綜合時間誤差和由此引入的綜合測量誤差.

表1 算法綜合誤差Table 1 Total error of the algorithm

從表1的實驗結(jié)果可以看出，在平均誤差方面，算法的綜合誤差隨著落體下落距離的增加呈現(xiàn)明顯的減小趨勢，但是均非常接近與零.而算法在累積、最大和最小綜合誤差方面，從2—18cm和6—18cm兩個下落段的計算結(jié)果中看，沒有明顯的改善.這從一個方面證明了，在干涉帶頻率比較低時，雖然算法誤差較大，但是由于過零點個數(shù)不多，算法的誤差沒有給結(jié)果帶來明顯的計算誤差.因此，可以利用該算法從落體下落2cm時開始取數(shù)計算，盡量增加有效的數(shù)據(jù)量，提高測量結(jié)果的可靠性.

3 實驗觀測結(jié)果與結(jié)論

2011年5—6月，筆者所屬研究小組使用基于本文介紹的軌跡重建算法的實驗樣機在北京國家地球觀象臺進行了試驗觀測.為研究落體次數(shù)對觀測精度的影響，在該臺的國家重力基準點（編號H03，基巖點）分別進行了320次下落和160次下落的兩次觀測.同時作為抗干擾觀測試驗，在該臺的普通辦公室設(shè)置重力點（編號H09，地板點）也進行了一次320次下落的測量.H03測點參考值采用FG5/232在2009年和2010年的觀測結(jié)果的平均值，H09測點的參考值采用兩套高精度相對重力儀（聯(lián)測精度10×10－8m/s2）與H03測點聯(lián)測得到.H09和H03兩點的參考重力值均根據(jù)梯度測量結(jié)果歸算到實驗樣機0.6m的有效高度.表2給出了這2個測點的3次測量結(jié)果與各測點的參考值之間的偏差與測量精度.

表2 2011年5月實驗樣機測量結(jié)果Table 2 Measurement results in the test of May 2011

結(jié)果表明，基于本文算法的實驗樣機精度在160次測量時即優(yōu)于10×10－8m/s2，系統(tǒng)偏差約71×10－8m/s2，H09測點兩次測量系統(tǒng)偏差固定.這一結(jié)果為下一步進行系統(tǒng)差的確認和研發(fā)更高精度的絕對重力儀提供了技術(shù)保障.

本文基于自主研制的數(shù)字化絕對重力儀的數(shù)字化干涉帶信號，提出了一種全新的數(shù)字處理算法，它不同于目前國際上流行的基于硬件的過零點觸發(fā)法和基于數(shù)字化干涉帶局部擬合算法.本文算法主要針對數(shù)字化的干涉帶信號，利用分段和對每段第一個過零點精確時刻的求取，提取落體下落過程中的時間位移坐標.然后通過最小二乘擬合求解重力加速度.可以說是基于硬件的過零觸發(fā)法的軟件替代.該設(shè)備和算法簡單，誤差小，采集可以實現(xiàn)模塊化，利于系統(tǒng)誤差的理論計算和固定，便于具有自主知識產(chǎn)權(quán)儀器的批量商業(yè)化生產(chǎn).

程力，張雅杰，蔡體菁.2007.重力輔助導(dǎo)航匹配區(qū)域選擇準則［J］.中國慣性技術(shù)學(xué)報，15（5）：559－－563.

郭有光，鐘斌，邊少峰.2003.地球重力場確定與重力場匹配導(dǎo)航［J］.海洋測繪，23（5）：61－－64.

劉棟勛，劉文泰.1984.重力預(yù)報地震動態(tài)［J］.四川地震，（1）：27－－31.

錢進，劉忠有，張小平，石春英，劉秀英，王悍平，蔡山，譚慧萍.2008.一種新型的碘穩(wěn)定633 nm He－Ne激光系統(tǒng)［J］.計量學(xué)報，29（1）：10－－13.

游澤霖，徐菊生.1984.精密絕對重力測量在地震研究中的應(yīng)用［J］.地殼形變與地震，4（4）：410－－414.

張善法，孟令順，杜曉娟，張鳳旭，馬國慶.2009.高精度重力測量在金礦采空區(qū)探測中的應(yīng)用研究［J］.地球物理學(xué)進展，24（2）：590－－595.

Bich W，D’Agostino G，Germak A，Pennecchi F.2007.Uncertainty propagation in a non－linear regression analysis：Application to a ballistic absolute gravimeter（IMGC－02）［C］∥Advanced Methods for Uncertainty Estimation in Measurement，2007 IEEE International Workshop：30－－34.

D’Agostino G，Germak A，Desogus S，Origlia C，Barbato G.2005.A method to estimate the time－position coordinates of a free－falling test－mass in absolute gravimetry［J］.Institute of Physics Publishing，42（4）：233－－238.

D’Agostino G，Desogus S，Germak A，Origlia C，Quagliotti D，Berrino G，Corrado G，d’Errico V Ricciardi G.2008.The new IMGC－02 transportable absolute gravimeter：Measurement apparatus and applications in geophysics and volcanology［J］.Annals of Geophysics，51（1）：39－－49.

de Angelis M，Bertoldi A，Cacciapuoti L，Giorgini A，Lamporesi G，Prevedelli M，Saccorotti G，Tino G M.2009.Precision gravimetry with atomic sensors［J］.Measurement Science and Technology，20（2）：1－－16.

Durando G，Desogus S，Mazzoleni F.2002.Accuracy assessment of analysis of a free－fall gravimeter data［C］∥Conference Digest 2002/Conference on Precision Electromagnetic Measurements.Ottawa，Ontario，Canada：86－－87.

Hinderer J，F(xiàn)lorsch N，Makinen J，Legros H，F(xiàn)aller J E.1991.On the calibration of a superconducting gravimeter using absolute gravity measurements［J］.Geophys J Int，106（2）：491－－497.

Rothleitner Ch，Svitlov S，Merimeche H，Hu H，Wang L J.2009.Development of new free－fall absolute gravimeters［J］.Metrologia，46（3）：283－－297.

Sakuma A.1984.Present status of the absolute measurement of gravitational acceleration［G］∥Taylor B N，Phillips W D eds.Precision Measurement and Fundamental ConstantsⅡ.Washington，D C：U S Government Printing Office：397－－404.

Simon D P，Williams，T F.Graham Jeffries B.2001.Absolute gravity measurements at UK tide gauges［J］.Geophys Res Lett，28（12）：2317－－2320.

Sugihara M.2003.Geothermal reservoir monitoring with an FG5 absolute gravimeter［C］∥Proceedings，Twenty－Eighth Workshop on Geothermal Reservoir Engineering.Stanford University，Stanford，California.

Vilibic I.1997.Global sea level rise？New techniques for the absolute sea level measurement［J］.Geofizika，14：119－－131.

Van Camp M，Camelbeeck T，Richard P.2003.The FG5 absolute gravimeter：Metrology and geophysics［J］.Physicalia Magazine，Journal of the Belgian Society，25（3）：161－－174.

Zumberge M A，Rinker R L，F(xiàn)aller J E.1982.A portable apparatus for absolute measurements of the earth’s gravity［J］.Metrologia，18（3）：145－－152.

A new type of algorithm for rebuilding the trace of free－fall body in absolute gravimeter development

Wu Qiong1），Teng Yuntian1）Huang Dalun2）Long Jianfeng1）
1）Institute of Geophysics，China Earthquake Administration，Beijing 100081，China
2）National Institute of Metrology P.R.China，Beijing 100013，China

10.3969/j.issn.0253－3782.2012.04.011

P315.61

A

中國地震局地球物理研究所基本科研業(yè)務(wù)費（DQJB01A02）資助.

2011－06－23收到初稿，2012－01－04決定采用修改稿.

e－mail：wuqiong＿cea＠hotmail.com

吳瓊，滕云田，黃大倫，龍劍鋒.2012.絕對重力儀研制中一種新的自由落體軌跡重建算法.地震學(xué)報，34（4）：549－－556.

Wu Qiong，Teng Yuntian，Huang Dalun，Long Jianfeng.2012.A new type of algorithm for rebuilding the trace of free－fall body in absolute gravimeter development.Acta Seismologica Sinica，34（4）：549－－556.