基于系數(shù)搜索的振動補償方法*

要佳敏 莊偉 馮金揚 王啟宇 趙陽 王少凱 吳書清 李天初

(中國計量科學(xué)研究院,北京 100029)

1 引言

絕對重力測量,即對重力加速度絕對值的精確測量,在輔助導(dǎo)航、資源勘探、地球物理、計量領(lǐng)域都有著廣泛的應(yīng)用[1-4].現(xiàn)有的高精度絕對重力儀普遍采用自由落體法,主要包括以角錐棱鏡為下落物體的激光干涉式絕對重力儀[5-8]和以冷原子團為下落物體的原子干涉式絕對重力儀[9-14],以下分別簡稱為光學(xué)重力儀和原子重力儀.前者通過激光干涉儀測量角錐棱鏡在真空中自由下落的運動軌跡,擬合求解重力加速度,后者根據(jù)原子團在自由下落的過程中經(jīng)過3 次拉曼激光脈沖作用形成的原子干涉條紋來計算重力加速度.由于二者的實際測量對象是下落物體或下落原子團相對于儀器內(nèi)一個參考鏡(角錐棱鏡或平面鏡)的運動加速度,因此參考鏡自身的加速度必然耦合到測量結(jié)果中.理論上,忽略參考鏡自身振動時,現(xiàn)有的絕對重力儀的測量精度可以達到微伽(μGal,1 μGal=1 ×10—8m·s—2)量級;但實際上,如果直接將參考鏡放置在地面上,地面振動引起的測值離散度將超過10 μGal,在復(fù)雜振動環(huán)境中甚至達到mGal 或10 mGal 量級.因此,利用測振儀器獲得參考鏡自身的運動加速度并在計算重力加速度時對其影響進行修正是提高絕對重力儀測量精度的有效手段之一,這種修正方法一般稱為振動補償.

具體而言,影響原子重力儀的地面振動噪聲主要為環(huán)境噪聲,與測量時間和測量點的地基、坐標(biāo)及環(huán)境有關(guān),可以根據(jù)頻率來粗略區(qū)分.首先是人類活動噪聲,一般大于1 Hz,如儀器測試期間附近的人員走動,車輛等人工振源產(chǎn)生的振動,以及建筑物的晃動等.其次是來自地球本身的噪聲,一般在0.1—1.0 Hz 之間,主要包括周期性的脈動和地震等地球內(nèi)部運動引起的振動,其中地脈動噪聲模型的加速度功率譜密度在頻率為0.2 Hz 和3.0 Hz的位置存在峰值[15].最后是氣壓波動等大氣運動噪聲,一般小于0.1 Hz,普遍變化緩慢且幅值較低,對絕對重力測量的影響可以忽略[16,17].由于絕對重力測量只依靠垂直方向的位移數(shù)據(jù)進行計算,因此對原子干涉條紋進行振動補償時,一般忽略地面振動的水平分量,僅需使用測振儀器輸出的垂直振動信號.

目前國際上較為成熟的原子重力儀是法國巴黎天文臺研制的CAG-01 型原子重力儀,曾參加過多次國際重力比對,采用的振動處理方法為被動式垂直隔振系統(tǒng)與振動補償?shù)慕Y(jié)合,其中振動補償?shù)木唧w原理為:數(shù)據(jù)采集卡同步記錄原子干涉過程中的地震計輸出的速度信號與原子干涉條紋,在控制軟件中利用無限沖激響應(yīng)(infinite impulse response,IIR)濾波器和非因果低通濾波器對地震計信號進行處理,以減小地震計傳遞函數(shù)幅頻特性不恒定、相頻特性存在非線性的影響,獲得更真實的參考鏡振動信號,最后利用處理后的信號得到由參考鏡振動引入的相位噪聲,修正原子干涉條紋.在該儀器所在的測量環(huán)境中,此振動補償方法本身可以使重力儀靈敏度提升為原來的3 倍[11].德國漢諾威大學(xué)提出利用數(shù)字式低通和高通濾波器對振動傳感器的輸出信號進行前處理,以得到更為精確的參考鏡振動引入的相位噪聲,從而實現(xiàn)振動補償.其實驗結(jié)果表明,在較安靜的環(huán)境中利用高精度商用地震計進行振動補償,或在復(fù)雜振動環(huán)境中使用精度較低的商用加速度計或其自主研制的新型光學(xué)慣性傳感器進行補償,均可以有效提高重力儀的靈敏度[18].在我國,浙江大學(xué)自制的原子重力儀同樣采用了振動補償方法,其主要特點為利用商用反演軟件修正地震計實測傳遞函數(shù)的影響以獲得更真實的參考鏡振動信號.在其實驗室環(huán)境中,該方法可以將重力儀的靈敏度提升至與使用被動隔振系統(tǒng)時相近的水平[19].國防科技大學(xué)對參考鏡垂向振動及水平偏轉(zhuǎn)引入的測量噪聲進行了詳細分析,利用加速度計測量原子重力儀工作時的參考鏡振動,同時對加速度計與參考鏡之間的傳遞函數(shù)進行了較為精確的標(biāo)定.實測結(jié)果表明在較復(fù)雜的振動環(huán)境中其振動補償方法可以實現(xiàn)顯著的振動補償效果[20].軍事科學(xué)院國防創(chuàng)新院提出的振動補償方法則采用具有高精度四通道移相探測器的邁克耳孫激光干涉儀替代傳統(tǒng)振動傳感器.該干涉儀可以直接測量出原子干涉條紋的相移,并結(jié)合信號整形、數(shù)字鑒相等方法確定拉曼激光脈沖作用前后原子干涉條紋的相位小數(shù)的變化,得到參考鏡振動引起的原子干涉相位偏差.去除該相位偏差,即可對重力儀測量到的躍遷概率信號進行修正[21].此外,中國計量科學(xué)研究院、中國科學(xué)院精密測量研究院、華中科技大學(xué)、中國科技大學(xué)、國防科技大學(xué)、美國斯坦福大學(xué)、德國洪堡大學(xué)等單位也研制了不同類型的隔振系統(tǒng)[12,13,22-26],以減小地面振動對原子重力儀的影響.

德國漢諾威大學(xué)的Richardson 等[18]和我國國防科技大學(xué)的研究人員[20]均指出,傳統(tǒng)振動傳感器與參考鏡由于安裝位置不同形成的機械傳遞函數(shù)可能是影響振動補償精度的重要因素;在原子重力儀進行正式測量前對該傳遞函數(shù)中的關(guān)鍵系數(shù)進行初始標(biāo)定有助于提高振動補償效果.這些分析為提升原子重力儀的振動補償精度提供了重要指導(dǎo),但目前尚無具有普適性的、能夠在任意測量時間及地點快速實現(xiàn)對傳遞函數(shù)進行實時修正的具體振動補償算法.而在光學(xué)重力儀方面,清華大學(xué)的研究人員給出了對重力測量結(jié)果進行實時修正的基于相同簡化模型的具體振動補償算法,并提供了該算法在不同地點對延時系數(shù)和增益系數(shù)進行實時計算的結(jié)果及重力測量的實測結(jié)果[27-29].他們采用的傳感器為高精度地震計,利用相關(guān)分析法或黃金分割法搜索出當(dāng)前環(huán)境下這2 個系數(shù)的具體值,并對下落物體軌跡信號進行修正.對同一套振動補償算法,輸入在中國計量科學(xué)研究院昌平院區(qū)的安靜環(huán)境中測得的重力數(shù)據(jù)和地震計數(shù)據(jù),算法的直接輸出結(jié)果表明振動補償后光學(xué)重力儀的測量分辨率可以提升5 倍;輸入在清華大學(xué)所處的復(fù)雜振動環(huán)境中的數(shù)據(jù),結(jié)果表明振動補償可將重力儀的測量分辨率提升16 倍,效果更為顯著[27].該算法最大的特點是較強的自適應(yīng)性,在不同時刻和不同地點的測量過程中可以直接使用,以保證真實的機械傳遞函數(shù)因環(huán)境變化而發(fā)生改變時對其當(dāng)前值進行實時估算,實現(xiàn)當(dāng)時當(dāng)?shù)刈詈玫恼駝友a償效果.受上述研究結(jié)果的啟發(fā),基于上述簡化模型,本文給出了一套振動補償算法的工作原理和具體流程,該振動補償算法可以應(yīng)用于原子重力儀，對傳遞函數(shù)進行實時標(biāo)定.我們通過仿真運算對算法的有效性進行了詳細驗證,最后利用實驗室自主研制的原子重力儀實際評估了該算法對原子干涉條紋擬合精度和重力測值離散度的補償效果,為提高原子重力儀在不同環(huán)境甚至是復(fù)雜振動環(huán)境中的測量精度提供幫助.

2 原 理

2.1 應(yīng)用于原子重力儀的振動補償原理

原子重力儀基于冷原子物質(zhì)波干涉原理,將激光脈沖作用到冷原子團,通過雙光子受激拉曼躍遷或多光子布拉格衍射過程來實現(xiàn)原子干涉過程[9].以基于受激拉曼躍遷的原子重力儀為例,冷原子團在真空腔中自由下落,真空腔下方的平面參考鏡將垂直向下的拉曼激光原路反射,使兩束對射的激光同時作用于原子團,在此過程中3 個連續(xù)的激光脈沖實現(xiàn)了原子波包的分束、反轉(zhuǎn)和合束,分別為第1 個π/2 脈沖、π 脈沖和第2 個π/2 脈沖,相鄰脈沖的時間間隔為T[30].合束后任一能態(tài)的原子均由2 條路徑上的原子疊加而成,通過熒光信號可以探測脈沖作用后原子團的躍遷概率P,即干涉條紋,滿足

式中:ΔФ為原子相位,與當(dāng)?shù)氐闹亓铀俣萭有關(guān);代表條紋對比度的系數(shù)A和B可以通過對干涉條紋進行余弦擬合得到.理想情況下,原子相位ΔФ等于其理論值ΔФth,即

式中,keff為有效波矢,α表示激光的掃頻速率.重力儀的控制軟件調(diào)控 (2) 式中的α線性增大,則原子相位ΔФ也將線性變化,使干涉條紋呈現(xiàn)為標(biāo)準(zhǔn)余弦信號的形式.此時P應(yīng)與由P擬合出的標(biāo)準(zhǔn)信號Pfit完全相同,有

但實際情況下,受到各種相位噪聲的影響,原子相位ΔФ等于其實際值ΔФm,有

式中Δφvib為參考鏡振動引入的相位噪聲,以下簡稱為振動相位;Δφothers為其他噪聲源引入的相位噪聲的總和,以下簡稱為其他相位噪聲.此時P由Pfit變?yōu)?/p>

當(dāng)振動相位在相位噪聲中占主導(dǎo)時,依據(jù)上述原理可以確定修正干涉條紋P的振動補償過程.第一步,對原始實測條紋P(ΔФth)進行余弦擬合,得到擬合出的理論條紋Pfit(ΔФth)和(1)式中的系數(shù)A和B.需要說明的是,此時可以將原子重力儀在正式測量前的初始測量過程中獲得的粗測值gth代入(2)式以得到擬合條紋Pfit.這是因為原子重力儀到達一個新的測量地點后一般都需要在啟動階段通過改變拉曼間隔T并掃描掃頻速率α來獲得正式測量所用的中心值αcenter,而該測點的粗測重力值即為gth=αcenter/keff.此后每一次或每一組躍遷概率數(shù)據(jù)都對應(yīng)真實重力值相對于該粗測值的偏差Δg,所以實際上重力儀最終輸出的高精度測值為g=gth+Δg.因此利用粗測值gth計算擬合條紋有利于提高計算效率,這也符合原子重力儀的一般測試流程[19,20].第二步,通過振動傳感器獲得一段時間t內(nèi)參考鏡的真實振動速度vm(t),進而根據(jù)

計算出真實的振動相位Δφvib,其中S(t)為重力儀的靈敏度函數(shù),T1和T3分別為第1 個π/2 脈沖和第2 個π/2 脈沖的作用時刻,均由重力儀控制軟件的已知參數(shù)決定.第三步,計算出僅受振動影響時的原子相位ΔФv和躍遷概率Pv,即

此時重新以ΔФv為橫坐標(biāo)、P為縱坐標(biāo)得到的數(shù)據(jù)對P(ΔФv)即為修正后的干涉條紋,至此補償過程結(jié)束.如圖1 所示,與藍色圓圈表示的原始條紋P(ΔФth)相比,紅色方塊表示的修正條紋P(ΔФv)的橫坐標(biāo)有所移動,理論上應(yīng)更接近灰色的標(biāo)準(zhǔn)余弦信號,即擬合出的理論條紋Pfit(ΔФth).也就是說,對修正后的條紋P(ΔФv)再次進行余弦擬合時的均方根誤差(root mean square error,RMSE)值應(yīng)小于對原始條紋P(ΔФth)進行擬合時的RMSE 值.

圖1 原子干涉條紋示意圖Fig.1.Schematic diagram of the fringe signal of atomic interferometer.

2.2 基于傳遞函數(shù)簡化模型的振動補償方法

基于傳遞函數(shù)簡化模型以實現(xiàn)上述振動補償過程的方法包括硬件和軟件兩方面.

硬件方面,由于目前尚無以原子重力儀的參考鏡為敏感質(zhì)量的地震計,因此一般而言只能用一臺獨立的地震計或加速度計放在參考鏡旁邊或下方測量其振動.然而,受機械結(jié)構(gòu)及地震計自身性能的限制,地震計輸出信號Us(t)并不能完全真實地反映參考鏡的振動速度vm(t),二者的關(guān)系如圖2所示,其中Ga表示從地面振動速度到參考鏡振動速度的傳遞函數(shù),Gb表示從地面振動速度到地震計敏感質(zhì)量振動速度的傳遞函數(shù),Gc表示從敏感質(zhì)量運動速度到地震計輸出電壓之間的傳遞函數(shù).由于地震計內(nèi)部的敏感質(zhì)量與參考鏡的實際位置必然存在水平和垂直差異,因此Ga與Gb不可能相等,且當(dāng)?shù)鼗牧细淖儠r,這兩個傳遞函數(shù)本身也可能發(fā)生變化.另一方面,地震計自身的傳遞函數(shù)Gc并非幅頻特性恒定、相頻特性為線性的理想傳遞函數(shù)[31],且有可能隨測量時間與環(huán)境的改變而發(fā)生變化(這是因為溫度、濕度、氣壓等環(huán)境參數(shù)會對地震計內(nèi)部的機電模塊參數(shù)產(chǎn)生影響,且這些參數(shù)可能隨時間發(fā)生漂移).綜上,參考鏡的真實振動速度與地震計輸出的電壓信號之間的傳遞函數(shù)H滿足

圖2 地震計輸出信號與參考鏡真實振動的關(guān)系Fig.2.Transfer function between the vibration of reference mirror and the output signal of seismometer.

因此對傳遞函數(shù)H的求解將直接影響振動補償?shù)男Ч?

雖然實際情況中Ga,Gb,Gc不易測量且可能改變,導(dǎo)致測量人員無法得到總傳遞函數(shù)H的精確值,但由于單次重力測量任務(wù)中測量環(huán)境不會有顯著變化,且測量時間一般不超過24 h,因此可以利用由一個延時系數(shù)τ和一個增益系數(shù)K組成的簡化模型來估算H,并認為這兩個系數(shù)的數(shù)值在當(dāng)前測量過程中保持不變[27].由此可得圖2 的簡化模型為

式中,Ks為地震計的標(biāo)稱靈敏度,由制造商提供.

理論上當(dāng)延時系數(shù)τ和增益系數(shù)K取最接近真實值時,利用地震計輸出電壓Us和(9)式計算出的參考鏡振動速度vm(t)也最接近真實值,對條紋的修正效果最好.因此振動補償可以選取對條紋進行擬合時的RMSE 值作為優(yōu)化目標(biāo),基本思路為:給定τ的取值范圍,遍歷搜索使RMSE 值達到最小值時對應(yīng)的最優(yōu)延時系數(shù)τopt;再給定K的取值范圍,在設(shè)定延時為τopt的情況下,再次遍歷搜索使RMSE 達到最小值時對應(yīng)的最優(yōu)增益系數(shù)Kopt.最后,在設(shè)定延時為τopt、增益為Kopt的情況下,得到修正后的條紋P(ΔФv)及其余弦擬合時的RMSE 值.

由此得到具體的算法流程如圖3 所示,相應(yīng)的處理程序即為振動補償方法的軟件部分.上述分析中使用的主要物理符號匯總在表1 中.

表1 主要物理符號描述Table 1.Description of main symbols.

圖3 基于系數(shù)搜索的振動補償算法流程Fig.3.The algorithm of the vibration correction method based on element searching.

2.3 振動補償前后干涉條紋擬合殘差的變化

上述振動補償算法的優(yōu)化目標(biāo)為擬合干涉條紋時的RMSE 值.該值具有簡單直觀的數(shù)學(xué)含義,但物理意義不明顯.類比基于相同簡化模型的振動補償方法在光學(xué)重力儀中的算法流程[27],可以用另一種更具有物理意義的變量作為振動補償效果的最終評價指標(biāo),即干涉條紋擬合殘差的變化.

運行振動補償算法前,先計算出原始條紋P(ΔФth)相較于理論條紋Pfit(ΔФth)的殘差Rm(ΔФth)的標(biāo)準(zhǔn)差σm,即

運行振動補償算法后,利用搜索出的最優(yōu)延時系數(shù)τopt和增益系數(shù)Kopt,根據(jù)(6)式,(7)式和(9)式得到僅受振動影響的干涉條紋Pv相較于理論條紋Pfit的殘差Rv.同時,記原始條紋P與僅受振動影響的條紋Pv的差值為Rc.顯然Rc也就是原始條紋殘差Rm與僅受振動影響的條紋殘差Rv的差值,相當(dāng)于去除振動影響后的條紋殘差.理論上振動補償效果越好,則振動相位Δφvib的計算越準(zhǔn)確,推算出的僅由振動引入的條紋殘差Rv也越接近原始條紋殘差Rm.此時Rc將趨近于僅由其他相位噪聲Δφothers引入的條紋殘差,其標(biāo)準(zhǔn)差σc也應(yīng)趨于最小值.因此可以以

相較原始條紋殘差Rm的標(biāo)準(zhǔn)差σm的衰減程度作為衡量振動補償效果的最終指標(biāo).

3 仿真驗證

上述振動補償算法的效果將在實驗室自主研發(fā)的NIM-AGRb-1 型原子重力儀[12]上進行考察.實際測量之前,可以利用美國 MathWorks 公司生產(chǎn)的 MATLAB 軟件驗證算法的有效性,具體的仿真流程為:1) 參考NIM-AGRb-1 型重力儀的控制參數(shù)和之前的實測數(shù)據(jù),給定合理的計算參數(shù)值,包括波矢keff、粗測重力值gth、靈敏度函數(shù)S(t)、拉曼激光的脈沖間隔T=80 ms、1 min 內(nèi)線性增大的掃頻速率序列α(時間間隔為2 s,相當(dāng)于原子團的下落間隔)、決定干涉條紋對比度的系數(shù)A和B等;2) 給定1 min 內(nèi)代表隨機振動的速度信號的幅值A(chǔ)rand和代表地脈動的速度信號的頻率序列fs及相應(yīng)的幅值序列As,具體取值參考地脈動噪聲功率譜密度的新高噪聲模型 (new high noise model,NHNM)[15],設(shè)定參考鏡真實振動速度信號vm(t)為這些信號的總和;3) 根據(jù)(6)式計算參考鏡振動引入的真實相位噪聲Δφvib,再根據(jù)(4)式和(5)式計算其他相位噪聲Δφothers為0 時的躍遷概率Psim,該概率Psim即代表原子重力儀輸出信號P的仿真值,數(shù)據(jù)量N=30,如圖4(a)所示;4)設(shè)定真實延時系數(shù)τset=5 ms,真實增益系數(shù)Kset=0.9,設(shè)定地震計靈敏度Ks為將在實測中使用的CMG-3 ESP 型地震計的靈敏度,根據(jù)(9)式反推出此時地震計敏感質(zhì)量的振動速度和地震計輸出電壓的仿真值vs(t)和Us(t),注意Us(t)即代表地震計在原子重力儀工作時同步輸出的信號,如圖4(b)所示;5)對上述仿真出的概率Psim和地震計電壓信號Us(t)運行振動補償算法.

圖4 (a) 仿真干涉條紋;(b) 仿真振動信號Fig.4.(a) The simulated fringe signal;(b) the simulated output of seismometer.

記以上仿真過程為第1 種情況,運算結(jié)果如圖5 和圖6 所示.可以看出,不考慮其他相位噪聲時,該算法可以精確地搜索出設(shè)定的延時系數(shù)τ=5 ms 和增益系數(shù)K=0.9(如圖中黑色虛線所示,下同),并將原始條紋Psim(ΔФth)調(diào)整為更接近理論條紋Pfit(ΔФth)的修正條紋Psim(ΔФv),其中典型的修正位置用橘色圓圈標(biāo)出.補償前對原始條紋Psim(ΔФth)進行余弦擬合時的RMSE 值為13.1×10—3,而補償后對修正條紋Psim(ΔФv)進行余弦擬合時的RMSE 值僅為0.02×10—3.同時該算法計算出的僅受振動影響的條紋相對于理論條紋的殘差Rv(ΔФth)和原始條紋相對于理論條紋的殘差Rm(ΔФth)完全相同.依據(jù)2.3 節(jié)提出的評價指標(biāo),此時原始條紋殘差Rm的標(biāo)準(zhǔn)差σm為12.39×10—3,而去除振動影響后的條紋殘差Rc的標(biāo)準(zhǔn)差σc僅為0.02×10—3,衰減幅度達到99.8%,驗證了該振動補償算法的有效性.

圖5 其他相位噪聲為0 時仿真得到的擬合RMSE 值與(a)延時系數(shù)和(b)增益系數(shù)的關(guān)系Fig.5.The RMSE value of fringe fitting as a function of (a) time delay and (b) proportional element,without considering other phase noises,calculated by simulation.

圖6 其他相位噪聲為0 時仿真得到的(a)修正前后的干涉條紋與理論擬合曲線(典型修正點用橘色圓圈標(biāo)出),(b)原始條紋的擬合殘差與計算出的僅由振動導(dǎo)致的條紋的擬合殘差Fig.6.(a) The simulated fringe signals before and after vibration correction with the fitted curve as the reference (with the typical corrected data marked by orange circles),and (b) the fitting residual of the simulated fringe signals,without considering other phase noises.

進一步,修改前述仿真流程的第3 步,設(shè)定其他相位噪聲Δφothers為±10 mrad 范圍內(nèi)的隨機值.這種設(shè)定可以更準(zhǔn)確地模擬實際測量,記該仿真過程為第2 種情況.此時雖然振動相位Δφvib的范圍約為±400 mrad,是其他相位噪聲Δφothers的40 倍,但仿真運算結(jié)果表明其他相位噪聲的干擾仍不可忽視.如圖7 和圖8 所示,此時算法搜索出的延時系數(shù)為τ=2.7 ms,與設(shè)定值的區(qū)別達到毫秒量級;增益系數(shù)為K=0.957,與設(shè)定值相比偏大6.3%.此時雖然該算法同樣可以將原始條紋Psim(ΔФth)調(diào)整為更接近理論條紋Pfit(ΔФth)的修正條紋Psim(ΔФv)(典型位置如橘色圓圈所示),但偶爾會因設(shè)定為隨機值的其他相位噪聲Δφothers而產(chǎn)生隨機運算誤差,如綠色方框所示的數(shù)據(jù)被修正到理論條紋曲線的另一側(cè)而不是曲線上;計算出的僅由振動引入的條紋殘差Rv(ΔФth)和原始條紋殘差Rm(ΔФth)大致相同,但也存在一定差異.最后,此時原始條紋殘差Rm的標(biāo)準(zhǔn)差σm為17.32×10—3,而去除振動影響后的條紋殘差Rc的標(biāo)準(zhǔn)差σc為7.82×10—3,衰減幅度為54.8%,明顯小于第1 種情況的衰減幅度.

圖7 考慮其他相位噪聲時仿真得到的擬合RMSE 值與(a)延時系數(shù)和(b)增益系數(shù)的關(guān)系Fig.7.The RMSE value of fringe fitting as a function of (a) time delay and (b) proportional element,with other phase noises considered,calculated by simulation.

圖8 考慮其他相位噪聲時仿真得到的(a)修正前后的干涉條紋與理論擬合曲線(典型修正點用橘色圓圈和綠色方框標(biāo)出),(b)原始殘差與計算出的僅由振動導(dǎo)致的條紋殘差Fig.8.(a) The simulated fringe signals before and after vibration correction with the fitted curve as the reference (with the typical corrected data marked by orange circles and green square),and (b) the fitting residual of the simulated fringe signals,with other phase noises considered.

上述結(jié)果說明其他相位噪聲將明顯限制振動補償算法的有效性.當(dāng)振動環(huán)境變差,即振動相位Δφvib在總相位噪聲中的權(quán)重增加時,算法的有效性將顯著提升.例如,仍取其他相位噪聲Δφothers為±10 mrad 范圍內(nèi)的隨機值,將參考鏡振動信號vm(t)中隨機振動的幅值A(chǔ)rand增大為原來的10 倍,并記該仿真過程為第3 種情況.此時擬合得到的延時系數(shù)為τ=4.1 ms,增益系數(shù)為K=0.919,比第2 種情況更接近設(shè)定值;原始條紋殘差Rm的標(biāo)準(zhǔn)差σm為29.28×10—3,去除振動影響后的條紋殘差Rc的標(biāo)準(zhǔn)差σc為6.78×10—3,衰減幅度為76.8%,高于第2 種情況.因此,理論上當(dāng)原子重力儀在復(fù)雜振動環(huán)境中測量時,該算法將發(fā)揮更為明顯的作用.

另一方面,增大數(shù)據(jù)量可以在小范圍內(nèi)提升振動補償算法對延時系數(shù)和增益系數(shù)的計算精度,但對條紋修正效果的提升水平有限.例如,設(shè)定其他相位噪聲Δφothers和振動相位Δφvib取值等同于第2 種情況,但將原子團下落的時間間隔減小為1 s.在仿真總時長仍為1 min 的情況下,此時躍遷概率信號的數(shù)據(jù)量擴大1 倍,達到N=60,記該仿真過程為第4 種情況.此時擬合得到的延時系數(shù)為τ=6.8 ms,增益系數(shù)為K=0.876,比第2 種情況更接近設(shè)定值;但原始條紋殘差Rm的標(biāo)準(zhǔn)差σm為14.80×10—3,去除振動影響后的條紋殘差Rc的標(biāo)準(zhǔn)差σc為6.56×10—3,衰減幅度為55.7%,和第2 種情況差別不大.

表2 匯總了以上所有運算結(jié)果.可以看出,補償前后對條紋進行余弦擬合時的RMSE 值的衰減幅度和條紋殘差標(biāo)準(zhǔn)差的衰減幅度及變化規(guī)律大致相同.如2.3 節(jié)所述,條紋殘差的變化更具有物理意義,因此在后續(xù)分析實測結(jié)果時不再專門提及RMSE 值的衰減幅度.

表2 仿真運算結(jié)果Table 2.Simulated results.

4 實驗結(jié)果

4.1 實驗設(shè)備

該振動補償算法的實際效果在中國計量科學(xué)研究院自主研制的NIM-AGRb-1 型原子重力儀上進行了驗證.NIM-AGRb-1 型重力儀具有兩種工作模式:其一為條紋掃描模式,在此期間拉曼激光波矢keff保持不變,掃頻速率α線性增加,可以得到原子干涉條紋;其二為重力測量模式,在此期間拉曼激光波矢keff不斷反向,即每一次作用于下落原子團時,激光脈沖的波矢keff均與上一次作用時數(shù)值相同、方向相反.此時,結(jié)合在原子鐘上運用的閉環(huán)鎖相技術(shù),可以實現(xiàn)對掃頻速率α的交叉鎖定,進而保證每一次原子團下落后都能得到重力測值g,而不必等待整個干涉條紋測量完畢后擬合出初相位進而計算出重力值,或依據(jù)基于不同拉曼間隔T得到的干涉條紋的交叉點位置來計算重力值.配合清華大學(xué)研制的被動式零長彈簧結(jié)構(gòu)的隔振系統(tǒng),該重力儀可實現(xiàn)44 μGal·Hz—1/2的靈敏度,合成不確定度約為5.2 μGal[12].

驗證該振動補償算法時,NIM-AGRb-1 型重力儀在第一種條紋掃描模式下工作,此時重力儀的輸出信號為實測躍遷概率P而非重力測值g,且最多掃描3 個整周期的干涉條紋.本實驗的現(xiàn)場情況如圖9(a)所示,測試地點為中國計量科學(xué)研究院昌平院區(qū)的重力實驗樓.實驗過程中參考鏡直接放置在重力儀探測腔下方的地面上,旁邊放置英國Guralp 公司生產(chǎn)的CMG-3 ESP 型地震計,其輸出信號的靈敏度為Ks=2000 V·m—1·s[31].探測腔中冷原子團下落期間的激光脈沖的時間間隔為T=80 ms,激光波長約為780 nm,相鄰原子團下落的時間間隔為2 s.受重力儀自身軟件的控制,采集躍遷概率P的數(shù)據(jù)采集卡的采樣頻率為1 Hz,無法用于對地震計輸出電壓的采樣,因此使用第二臺美國National Instruments 公司生產(chǎn)的PXIe-6361 數(shù)據(jù)采集卡[32]來記錄地震計的輸出信號Us(t),采樣頻率為500 kS/s,如圖9(b)所示.該采集卡還同步采集了重力儀的某一路觸發(fā)信號Utrg(t).重力儀的控制軟件控制采集該信號總比原子團下落時的第1 個π/2 脈沖提前40 ms,從而確保兩臺采集卡分別采樣的數(shù)據(jù)可以在時間上互相匹配.

圖9 實驗現(xiàn)場圖 (a) NIM-AGRb-1 型原子重力儀;(b)地震計信號及觸發(fā)信號采集系統(tǒng)Fig.9.Photograph of (a) NIM-AGRb-1 atomic-interferometry absolute gravimeter,and (b) the data acquisition device recording the voltage signal of seismometer and trigger.

該測試地點位于地廣人稀的郊區(qū),周圍無交通、建筑和大型機械振源,且具有隔振地基,屬于非常安靜的測試環(huán)境,其地面振動加速度的功率譜密度如圖10 所示.由于NIM-AGRb-1 型重力儀體積較大,尚無法搬運至其他地點,因此本次實驗暫且僅在該安靜環(huán)境中進行.之后該算法的振動補償效果也將在振動條件更為復(fù)雜的環(huán)境中驗證.

圖10 測試地點的地面振動加速度功率譜密度Fig.10.The power spectrum density of the ground vibration acceleration at the measurement site.

4.2 實驗結(jié)果

本次實驗首先采集了10 組數(shù)據(jù),每組數(shù)據(jù)包含30 次原子團下落過程,其他控制參數(shù)如前所述.對每組數(shù)據(jù)分別運行振動補償算法,結(jié)果如表3 所示.

表3 10 組數(shù)據(jù)的實際測量結(jié)果Table 3.Measurement results obtained from 10 data sets.

可以看出該振動補償算法在安靜環(huán)境中已具有良好的條紋修正效果,去除振動影響后的條紋殘差Rc的標(biāo)準(zhǔn)差σc相比原始條紋殘差Rm的標(biāo)準(zhǔn)差σm的衰減幅度的均值達到44.4%,最大可達58.2%,與表1 中第2 種仿真情況的運算結(jié)果基本一致;相比對原始條紋進行余弦擬合時的RMSE值,擬合修正條紋時的RMSE 值也顯著減小.另一方面,各組數(shù)據(jù)分別計算出的最優(yōu)增益系數(shù)Kopt的標(biāo)準(zhǔn)差為0.115(均值為1.119),計算誤差較小;而最優(yōu)延時系數(shù)τopt的標(biāo)準(zhǔn)差為2.47 ms(均值為3.56 ms),計算誤差較大.前述表2 中第2 種情況的仿真結(jié)果已說明,即使Δφothers比Δφvib小1 個量級,τopt的計算誤差也將達到毫秒量級.因此最有可能導(dǎo)致實測結(jié)果中τopt的標(biāo)準(zhǔn)差達到毫秒量級的原因仍是其他相位噪聲Δφothers.這再次說明在后續(xù)工作中有必要在復(fù)雜振動環(huán)境中驗證該算法的振動補償效果,即考察振動相位Δφothers在總相位噪聲中所占的權(quán)重更大時,σc相比σm的衰減幅度以及延時系數(shù)和增益系數(shù)的計算精度能否如仿真預(yù)測一樣提升.

可以以第1 組數(shù)據(jù)的計算結(jié)果為例觀察振動補償算法的工作過程.如圖11 所示,該算法可以有效搜索出擬合條紋的RMSE 達到最小值時對應(yīng)的唯一最優(yōu)延時系數(shù)τopt和最優(yōu)增益系數(shù)Kopt.從圖12(a)可以看出,相比原始條紋P(ΔФth),對橫坐標(biāo)進行振動補償后的修正條紋P(ΔФv)更接近理論條紋Pfit(ΔФth)(典型位置如橘色圓圈所示).補償前對原始條紋P(ΔФth)進行余弦擬合時的RMSE 值為11.2 × 10—3,而補償后對修正條紋Psim(ΔФv)進行余弦擬合時的RMSE 值減小為4.8 × 10—3.原始條紋P相對于理論條紋的殘差Rm(ΔФth)和從地震計實測信號推算出的僅受振動影響時的條紋Pv相對于理論條紋的殘差Rv(ΔФth)也基本吻合,如圖12(b)所示.

圖11 第1 組實測數(shù)據(jù)的擬合RMSE 值與(a)延時系數(shù)和(b)增益系數(shù)的關(guān)系Fig.11.The RMSE value of fringe fitting as a function of (a) time delay and (b) proportional element,calculated from the 1st dataset of measurement.

圖12 第1 組實測數(shù)據(jù)(a)修正前后的干涉條紋與理論擬合曲線(典型修正點用橘色圓圈標(biāo)出),(b)原始殘差與計算出的僅由振動導(dǎo)致的條紋殘差Fig.12.(a) The fringe signals before and after vibration correction with the fitted curve as the reference (with the typical corrected data marked by orange circles),and (b) the fitting residual of the fringe signals,calculated from the 1st dataset of measurement.

實際上還可以采用基于黃金分割原理的二維搜索方法重新編寫振動補償算法流程,以實現(xiàn)對最優(yōu)延時系數(shù)τopt和最優(yōu)增益系數(shù)Kopt的全局搜索.計算結(jié)果表明,此時振動補償算法的效果與表3 所示的分別搜索τopt和Kopt時的效果相當(dāng).這也再次證明(9)式所示的傳遞函數(shù)簡化模型具有合理性.

基于表3 所列的10 組干涉條紋及后續(xù)測量得到的50 組干涉條紋進行計算的結(jié)果如圖13 所示,其中圖13(a)為這60 組條紋各自修正前的殘差標(biāo)準(zhǔn)差σm以及修正后的殘差標(biāo)準(zhǔn)差σc,圖13(b)為由此得到的σc相對于σm的衰減比.可以看出每組條紋殘差標(biāo)準(zhǔn)差的衰減比的均值仍為45%左右,且分布特征具有隨機性,導(dǎo)致該隨機性分布最可能的原因仍是原始數(shù)據(jù)中其他相位噪聲造成的隨機干擾.

圖13 修正前后60 組實測干涉條紋的(a)殘差標(biāo)準(zhǔn)差和(b)殘差標(biāo)準(zhǔn)差的衰減比Fig.13.(a) The standard deviations of the fitting residuals of 60 sets of fringe signals with and without the vibration correction,and (b) the reducing factor between them.

如前所述,進行本實驗時NIM-AGRb-1 型重力儀工作在條紋掃描模式,重力儀本身不輸出依據(jù)干涉條紋得到的重力測值g.此時,為了分析該振動補償算法對重力測值的影響,可以參考浙江大學(xué)提出的方法[19]計算重力測值的離散度:1)對每一組條紋進行余弦擬合時,進一步計算擬合曲線的初相位φinitial(n)(n為條紋序號,n=1,2,···);2)去除φinitial(n)的均值,得到初相位的偏差序列Δφinitial(n);3)根據(jù)(4)式,可得重力偏差序列

也就是說,基于修正前后的干涉條紋P(ΔФth)和P(ΔФvib),可以依據(jù)(12)式得到振動補償前后的重力偏差序列Δgm(n)和Δgc(n).由于本實驗所在的測量環(huán)境無固定振源,僅有與測量過程在時間上不相關(guān)的地脈動、輕微的人員走動或?qū)嶒灅峭獾能囕v移動等隨機振源,理論上參考鏡振動只影響重力測值的離散度,不會產(chǎn)生恒定的系統(tǒng)誤差.因此補償前后重力偏差序列Δgm(n)和Δgc(n)的對比可以合理地說明該振動補償算法的有效性.

根據(jù)上述分析,由修正前后的條紋數(shù)據(jù)得到的修正前后的重力偏差如圖14 所示.對全部60 組數(shù)據(jù),振動補償后重力測值的離散度相比補償前降低了一半左右;由修正前的實測條紋P(ΔФth)得到的重力偏差Δgm的標(biāo)準(zhǔn)差為σgm=51.5 μGal,而由修正后的條紋P(ΔФvib)得到的重力偏差Δgc的標(biāo)準(zhǔn)差σgc=22.0 μGal.σgc相對于σgm的衰減幅度達到57.3%,說明通過該振動補償算法提高干涉條紋的擬合精度可以顯著提升重力測量的精度.

圖14 由60 組實測干涉條紋得到的修正前后的重力偏差Fig.14.The variation of the measured g value deduced from 60 sets of fringe signals with and without the vibration correction.

目前NIM-AGRb-1 型原子重力儀作為計量基準(zhǔn)裝置已被設(shè)定為長期工作在第二種重力測量模式,即掃頻速率α處于鎖定狀態(tài),無法直接采用本文提出的振動補償算法進行測量.因此本文仍以考察該算法對原子干涉條紋殘差的衰減效果為主.目前實驗室正在研制采用條紋掃描模式來長期工作的可移動式原子重力儀,該振動補償算法也將配合此重力儀在不同的測試環(huán)境中進行長期測量,驗證其對包括阿倫方差、重力測量靈敏度在內(nèi)的關(guān)鍵重力參數(shù)的提升效果.

5 討論

在前述仿真和實驗過程中,振動補償算法查找最優(yōu)延時系數(shù)時的優(yōu)化目標(biāo)是擬合條紋時的RMSE 值,即取RMSE 達到其最小值時對應(yīng)的延時系數(shù)為最優(yōu)延時系數(shù)τopt.另一方面,根據(jù)理論分析及圖6(b)、圖8(b)和圖11(b)所示的仿真、實測結(jié)果,理論上振動補償效果好時,從地震計信號推算出的僅受振動影響時的條紋相對于理論條紋的殘差Rv(ΔФth)非常接近原始條紋相對于理論條紋的殘差Rm(ΔФth).因此,參考基于相同簡化模型的振動補償方法在光學(xué)重力儀中的算法流程[27],實際上還可以嘗試以Rv和Rm的相關(guān)系數(shù)CR作為優(yōu)化目標(biāo)搜索最優(yōu)延時系數(shù)τopt.

相關(guān)系數(shù)表示2 個列向量之間線性相關(guān)的程度,取值范圍為—1 至1,絕對值越大表示相關(guān)性越強[33].因此,相關(guān)系數(shù)CR的絕對值越大,說明Rv(ΔФth)和Rm(ΔФth)越接近,即推算出的振動相位Δφvib越準(zhǔn)確,則修正后的條紋越接近擬合曲線,其RMSE值應(yīng)趨于最小值.需要說明的是,計算相關(guān)系數(shù)CR時需要對Rc和Rv進行歸一化處理,因此無法以其為優(yōu)化目標(biāo)來搜索最優(yōu)增益系數(shù)Kopt.

以CR達到最大值時對應(yīng)的延時系數(shù)為最優(yōu)延時系數(shù)τopt,相應(yīng)的仿真和實測數(shù)據(jù)的計算結(jié)果分別如表4 和表5 所示.總體而言,采用兩種優(yōu)化目標(biāo)時分別計算出的最優(yōu)延時系數(shù)τopt非常接近.如圖15 所示,RMSE 值隨延時系數(shù)的變化曲線與相關(guān)系數(shù)CR隨延時系數(shù)的變化曲線具有幾乎完全相同的變化規(guī)律,但變化方向相反,符合理論預(yù)期.個別情況下采用兩種優(yōu)化目標(biāo)時分別計算出的最優(yōu)延時系數(shù)τopt存在毫秒量級的差異,例如表4中第4 種情況和表4 中第2 組數(shù)據(jù),其中前者對延時系數(shù)的遍歷搜索結(jié)果如圖15(d)所示,RMSE 值達到最小值時的延時系數(shù)和相關(guān)系數(shù)CR達到最大值時的延時系數(shù)略有區(qū)別(分別用黑色和粉色點劃線標(biāo)出).導(dǎo)致這種情況最可能的原因仍是原始數(shù)據(jù)中其他相位噪聲造成的隨機干擾.

表4 以相關(guān)系數(shù)CR 為優(yōu)化目標(biāo)時的仿真運算結(jié)果Table 4.Simulated results with the correlation value CR as the objective of the algorithm.

圖15 擬合RMSE 值和相關(guān)系數(shù)CR 與延時系數(shù)的關(guān)系 (a) 仿真時的第1 種情況;(b) 仿真時的第2 種情況;(c) 實測時的第1 組數(shù)據(jù);(d) 仿真時的第4 種情況Fig.15.The RMSE value of fringe fitting and the correlation value CR between the fitting residual of fringes before and after the vibration correction as a function of time delay of (a) 1st simulated dataset,(b) 2nd simulated dataset,(c) 1st measured dataset and(d) 4th simulated dataset.

表5 以相關(guān)系數(shù)CR 為優(yōu)化目標(biāo)時的10 組數(shù)據(jù)的實際測量結(jié)果Table 5.Measurement results obtained from 10 data sets with the correlation value CR as the objective of the algorithm.

對每一組數(shù)據(jù),依據(jù)振動補償算法計算出最優(yōu)延時系數(shù)τopt和增益系數(shù)Kopt后,可以得到最終確定的僅受振動影響時的條紋相對于理論條紋的殘差Rv,進而得到最終確定的相關(guān)系數(shù)CR.圖16展示了60 組數(shù)據(jù)中每一組數(shù)據(jù)的相關(guān)系數(shù)CR以及修正后的條紋殘差標(biāo)準(zhǔn)差σc相對于修正前的條紋殘差標(biāo)準(zhǔn)差σm的衰減比的對應(yīng)情況.可以看出,相關(guān)系數(shù)CR越大,σc相對于σm的衰減比也越大,即振動補償?shù)男Ч胶?二者的變化趨勢完全吻合,與理論預(yù)期相符,再次驗證了該振動補償算法原理的合理性.

圖16 最終確定的Rv 和Rm 的相關(guān)系數(shù)CR 與修正前后條紋殘差標(biāo)準(zhǔn)差的衰減比Fig.16.The correlation value CR between the fitting residual of fringes before and after the vibration correction,and the corresponding reducing factor of the standard deviation of the fitting residuals of fringe signals after the vibration correction.

還可以直接以2.3 節(jié)提出的條紋殘差標(biāo)準(zhǔn)差這一指標(biāo)作為振動補償?shù)膬?yōu)化目標(biāo),即取(11)式所示的標(biāo)準(zhǔn)差σc達到最小值時的延時系數(shù)和增益系數(shù)為最優(yōu)值.相應(yīng)的仿真及實測數(shù)據(jù)的計算結(jié)果表明,在絕大多數(shù)情況下以RMSE 值或σc為優(yōu)化目標(biāo)時搜索出的τopt和Kopt均完全相同.個別情況下搜索出的延時系數(shù)存在毫秒量級的差異,增益系數(shù)也可能存在微小差異,如表6 所示.導(dǎo)致這種結(jié)果最可能的原因仍是原始數(shù)據(jù)中其他相位噪聲造成的隨機干擾.

表6 以殘差標(biāo)準(zhǔn)差為優(yōu)化目標(biāo)時的仿真及實測數(shù)據(jù)運算結(jié)果Table 6.Simulated and measurement results with the standard deviation of the fitting residuals as the objective of the algorithm.

綜上,當(dāng)原子重力儀在振動條件更復(fù)雜的環(huán)境中測量時,即振動相位Δφvib的影響遠大于其他相位噪聲Δφothers的影響時,該振動補償算法的計算

精度和對條紋的修正效果有望得到進一步提升.在后續(xù)工作中,該振動補償算法將配合實驗室正在研制的可移動式原子重力儀在復(fù)雜環(huán)境中進行測試.

總體而言,本文在第2 節(jié)定性分析了原子重力儀所用參考鏡與振動傳感器之間的傳遞函數(shù)簡化模型中的延時系數(shù)τ和增益系數(shù)K發(fā)生變化的原因,描述了對這兩項系數(shù)進行搜索求解的具體振動補償算法流程;在第3 節(jié)采用不同的仿真輸入對算法的有效性進行了詳細的定量論證;在第4 節(jié)利用重力儀在非常安靜的環(huán)境中的實測數(shù)據(jù)驗證了該算法對原子干涉條紋以及重力測值的補償效果.該算法最大的優(yōu)點是自適應(yīng)性較強:相比現(xiàn)有的其他振動補償研究成果,它具有較強的技術(shù)意義,可以直接應(yīng)用于任何能夠掃描出干涉條紋的原子重力儀在任意時間和地點的測量結(jié)果,以實現(xiàn)對傳遞函數(shù)中的兩項關(guān)鍵系數(shù)的實時標(biāo)定,實現(xiàn)當(dāng)前最好的振動補償效果.另一方面,受限于NIM-AGRb-1型原子重力儀的體積、工作模式和測試任務(wù),該算法的有效性還未能在較為復(fù)雜的振動環(huán)境中進行驗證;而目前也尚無現(xiàn)有的其他振動補償方法在與本次實驗相同的安靜環(huán)境中應(yīng)用于與該重力儀精度水平相當(dāng)?shù)脑又亓x上的實驗數(shù)據(jù).從該算法與其他振動補償方法在不同測量環(huán)境及具有不同噪聲水平的原子重力儀上的實驗結(jié)果來看,該算法尚未展現(xiàn)出顯著優(yōu)于其他方法的補償效果,有待更進一步的測試驗證,其工作原理和具體流程也將在后續(xù)工作中不斷優(yōu)化.

6 結(jié)論

應(yīng)用于原子干涉式絕對重力儀的振動補償方法的主要作用是修正參考鏡振動引入的相位噪聲,從而提高干涉條紋的擬合精度,最終提高重力加速度g的測量精度.現(xiàn)有的振動補償方法在確定的靜態(tài)或動態(tài)環(huán)境中可以實現(xiàn)較高的測量精度,但大部分公開成果沒有介紹對于因測量環(huán)境改變而產(chǎn)生的原子重力儀中某些部件的機械傳遞函數(shù)的變化進行實時修正的具體應(yīng)對流程.從長期穩(wěn)定性的角度來看,這可能在一定程度上限制重力儀在不同環(huán)境尤其是復(fù)雜振動環(huán)境的測量精度.本文介紹了一種應(yīng)用于原子重力儀的基于搜索傳遞函數(shù)簡化模型系數(shù)的振動補償方法,著重描述了其具體的算法流程,并利用仿真及實測結(jié)果對其合理性和有效性進行了詳細論證.具體而言,該方法將參考鏡的振動信號與測量該振動的傳感器的實際輸出信號之間的真實傳遞函數(shù)簡化為一個延時系數(shù)和一個增益系數(shù)的結(jié)合,通過計算軟件對重力儀的干涉條紋和傳感器的輸出信號進行分析,搜索出當(dāng)前測量環(huán)境下這兩個系數(shù)的最優(yōu)估算值,由此可以更加精確地推算出當(dāng)前環(huán)境下參考鏡的真實振動信號.該方法具有自適應(yīng)性,可以在測量環(huán)境改變時對結(jié)構(gòu)的真實傳遞函數(shù)進行合理估算,以實現(xiàn)在當(dāng)時當(dāng)?shù)貙υ痈缮鏃l紋最好的修正效果.實測結(jié)果表明,在極為安靜的環(huán)境中,本文提出的振動補償算法可以將原子重力儀干涉條紋的余弦擬合殘差的標(biāo)準(zhǔn)差大幅衰減,最大衰減幅度達到58%.利用修正后的干涉條紋計算出的重力值的標(biāo)準(zhǔn)差相比修正前衰減57%.根據(jù)理論分析和仿真運算,在復(fù)雜振動環(huán)境中進行重力測量時,該算法將表現(xiàn)出更好的干涉條紋修正效果.因此,該振動補償算法有望直接應(yīng)用在不同時刻和不同地點對原子重力儀測量結(jié)果的實時修正,提高原子重力儀在不同環(huán)境尤其是復(fù)雜振動環(huán)境中的測量精度,其補償效果也將在之后的工作中進行更充分的驗證及優(yōu)化.