基于擴(kuò)展卡爾曼濾波算法的船載絕對(duì)重力測(cè)量數(shù)據(jù)處理*

朱棟 徐晗 周寅 吳彬 程冰 王凱楠 陳佩軍 高世騰 翁堪興 王河林 彭樹(shù)萍 喬中坤 王肖隆 林強(qiáng)

(浙江工業(yè)大學(xué)理學(xué)院，浙江省量子精密測(cè)量重點(diǎn)實(shí)驗(yàn)室，杭州 310023)

1 引言

基于冷原子干涉儀的精密測(cè)量技術(shù)近些年得到了快速發(fā)展.冷原子重力儀[1-6]逐漸進(jìn)入實(shí)用化階段，儀器的精度、靈敏度、穩(wěn)定性、尺寸、質(zhì)量、功耗、環(huán)境適應(yīng)性等性能都有了較大提升[2，4，7-10].這使得外場(chǎng)絕對(duì)重力測(cè)量得以快速發(fā)展[5，8，11-15].基于冷原子重力儀的絕對(duì)重力動(dòng)態(tài)測(cè)量實(shí)驗(yàn)也開(kāi)始出現(xiàn).國(guó)內(nèi)外已經(jīng)開(kāi)展了船載系泊絕對(duì)重力測(cè)量[16]、航空絕對(duì)重力測(cè)量[17]、海洋絕對(duì)重力測(cè)量[18]等實(shí)驗(yàn).傳統(tǒng)的重力場(chǎng)動(dòng)態(tài)測(cè)量均為相對(duì)測(cè)量，長(zhǎng)時(shí)間測(cè)量中存在零點(diǎn)漂移問(wèn)題，在動(dòng)態(tài)環(huán)境中基于原子重力儀的絕對(duì)重力測(cè)量可為其提供校準(zhǔn)，因此，絕對(duì)重力動(dòng)態(tài)測(cè)量具有重要的應(yīng)用價(jià)值.

在動(dòng)態(tài)條件下，測(cè)量噪聲的抑制對(duì)測(cè)量性能的提升至關(guān)重要，而動(dòng)態(tài)絕對(duì)重力數(shù)據(jù)的處理是其中的一個(gè)關(guān)鍵因素.測(cè)量噪聲的抑制除了在硬件方面提升原子重力儀整體的穩(wěn)定性和抗干擾能力外，還可以從算法角度進(jìn)行處理.原子重力儀通常的數(shù)據(jù)分析算法是通過(guò)最小二乘法對(duì)每條干涉條紋數(shù)據(jù)進(jìn)行余弦函數(shù)擬合[2，19]，使得每組條紋的估計(jì)值和觀測(cè)值之間的方差最小，從而獲得重力值的估計(jì).在實(shí)驗(yàn)室條件下，上述方法比較適用，但是對(duì)于海上的動(dòng)態(tài)測(cè)量環(huán)境而言，由于測(cè)量噪聲的增大，這種方法的適用性與實(shí)時(shí)性會(huì)變差.利用適用于船載測(cè)量系統(tǒng)的非線性濾波算法可降低測(cè)量數(shù)據(jù)的噪聲水平，從而獲得絕對(duì)重力值的估計(jì).目前，法國(guó)開(kāi)展的船載動(dòng)態(tài)絕對(duì)重力測(cè)量利用其自研的組合算法、原子重力儀與經(jīng)典加速度計(jì)之間的組合協(xié)議以及貝塞爾低通濾波進(jìn)行數(shù)據(jù)處理[18]，該數(shù)據(jù)處理方法雖然有效，但是在實(shí)時(shí)性上表現(xiàn)不強(qiáng).

本文提出基于擴(kuò)展卡爾曼濾波(extended Kalman filter，EKF)算法的數(shù)據(jù)處理方法對(duì)船載絕對(duì)重力測(cè)量數(shù)據(jù)進(jìn)行處理以提升實(shí)時(shí)性，并抑制測(cè)量噪聲.EKF 是一種基于線性卡爾曼濾波(最優(yōu)濾波器)的非線性時(shí)域?yàn)V波方法，適用于非線性系統(tǒng)的狀態(tài)估計(jì)(估計(jì)是通過(guò)系統(tǒng)的觀測(cè)值對(duì)系統(tǒng)的狀態(tài)量進(jìn)行定量推斷的過(guò)程)與追蹤[20-22]，其利用系統(tǒng)的狀態(tài)參量、狀態(tài)模型和觀測(cè)模型來(lái)描述系統(tǒng)隨時(shí)間的演變.在原子干涉式慣性導(dǎo)航傳感器[23，24]、原子傳感器的信號(hào)追蹤[25]、導(dǎo)航兼容的量子加速度計(jì)[26]等與量子傳感器相關(guān)的領(lǐng)域都有應(yīng)用.與最小二乘法的余弦擬合方法相比，EKF 算法不僅以遞推方式使方差最小，并且將系統(tǒng)狀態(tài)量的方差、均值等統(tǒng)計(jì)特性作為其參數(shù)，通過(guò)局部的線性化解析近似獲得系統(tǒng)狀態(tài)量的估計(jì).不僅可以有效追蹤船載絕對(duì)重力測(cè)量系統(tǒng)狀態(tài)量的變化，同時(shí)還可以實(shí)現(xiàn)重力測(cè)量數(shù)據(jù)的噪聲抑制.

本文根據(jù)船載冷原子干涉絕對(duì)重力測(cè)量系統(tǒng)的物理模型及其特點(diǎn)，引入了適用于該系統(tǒng)的EKF 算法，并對(duì)測(cè)量數(shù)據(jù)進(jìn)行濾波迭代處理，抑制了部分測(cè)量噪聲，獲得了重力值的估計(jì).航速小于2.1 km/h 條件下的絕對(duì)重力測(cè)量靈敏度得到了有效提升，并將處理得到的絕對(duì)重力數(shù)據(jù)與高精度地球重力模型XGM2019 計(jì)算的數(shù)據(jù)進(jìn)行了比對(duì)，發(fā)現(xiàn)兩者結(jié)果具備較好的一致性.本文為船載冷原子干涉式絕對(duì)重力測(cè)量系統(tǒng)在動(dòng)態(tài)測(cè)量條件下重力值的估計(jì)和重力測(cè)量噪聲的抑制提供了一種新的處理方法.

2 基于冷原子干涉儀的船載絕對(duì)重力動(dòng)態(tài)測(cè)量及數(shù)據(jù)處理原理

原子重力儀利用雙光子的受激拉曼躍遷，使銣87 原子在重力場(chǎng)下發(fā)生物質(zhì)波干涉，最后從干涉信號(hào)中提取出絕對(duì)重力值[2，5，14，27].船在海上航行時(shí)，原子重力儀的傾斜會(huì)隨船的運(yùn)動(dòng)不斷改變，需要利用慣性穩(wěn)定平臺(tái)實(shí)時(shí)維持重力傳感器在垂線方向.此外，海浪和船發(fā)動(dòng)機(jī)產(chǎn)生的振動(dòng)會(huì)傳遞到原子重力儀，其產(chǎn)生的額外相位噪聲會(huì)進(jìn)入到原子干涉儀的總相位中，導(dǎo)致無(wú)法得到有效的干涉條紋相位.根據(jù)相關(guān)文獻(xiàn)[16，18]，可以將振動(dòng)轉(zhuǎn)化為原子干涉時(shí)的相位干擾φvib，從而進(jìn)行修正.振動(dòng)引起的相位干擾經(jīng)過(guò)修正后，并不能完全消除其影響，還存在一部分的殘余相位噪聲.因此，原子干涉條紋的相位噪聲還包含其他一些難以定量分析的隨機(jī)噪聲.EKF 算法對(duì)船載冷原子干涉式絕對(duì)重力測(cè)量系統(tǒng)的狀態(tài)估計(jì)和噪聲處理是一種優(yōu)選的方法，因此，將其納入船載絕對(duì)重力測(cè)量系統(tǒng)的框架之內(nèi).

船載冷原子干涉式絕對(duì)重力測(cè)量系統(tǒng)是非線性系統(tǒng)，其大部分噪聲是高斯型噪聲，將EKF 算法作為該系統(tǒng)的估計(jì)器是一個(gè)不錯(cuò)的選擇.在狀態(tài)參量噪聲特性的合理統(tǒng)計(jì)評(píng)估下，可以有效追蹤和估計(jì)原子干涉條紋的相關(guān)參數(shù)，增大絕對(duì)重力估計(jì)數(shù)據(jù)隨時(shí)間和空間變化的實(shí)時(shí)性.通過(guò)狀態(tài)量的不確定性和測(cè)量噪聲之間的權(quán)衡，可以抑制重力測(cè)量噪聲.此外，EKF 算法以迭代遞推的方式進(jìn)行船載絕對(duì)重力測(cè)量系統(tǒng)狀態(tài)的濾波估計(jì)，需要的內(nèi)存較小、計(jì)算效率較高.

船載絕對(duì)重力動(dòng)態(tài)測(cè)量的原理如圖1 所示.其中，銣原子在真空腔中經(jīng)過(guò)冷原子團(tuán)的制備、原子選態(tài)、原子干涉等過(guò)程后，利用歸一化探測(cè)獲得原子布居數(shù)據(jù).通過(guò)線性?huà)呙枥獾念l率啁啾率獲取原子干涉條紋.經(jīng)典加速計(jì)采集的振動(dòng)信號(hào)與原子干涉儀的響應(yīng)函數(shù)進(jìn)行卷積獲得振動(dòng)相位干擾φvib.原子干涉條紋數(shù)據(jù)、拉曼光的頻率啁啾率和振動(dòng)的相位干擾被輸入到EKF 濾波算法以進(jìn)行重力值的估計(jì)，最后獲得噪聲水平更低的重力測(cè)量值.

圖1 船載冷原子干涉式絕對(duì)重力動(dòng)態(tài)測(cè)量原理Fig.1.The principle of absolute gravity dynamic measurement based on cold atom interferometer on ship.

根據(jù)船載絕對(duì)重力測(cè)量系統(tǒng)的物理模型可以建立該系統(tǒng)的觀測(cè)方程和狀態(tài)方程，從而對(duì)系統(tǒng)狀態(tài)的演化進(jìn)行描述.EKF 算法通過(guò)系統(tǒng)狀態(tài)的轉(zhuǎn)移和觀測(cè)過(guò)程的一階線性化解析近似，對(duì)當(dāng)前時(shí)刻的原子布居數(shù)和系統(tǒng)狀態(tài)量進(jìn)行預(yù)測(cè);同時(shí)，使用當(dāng)前時(shí)刻的觀測(cè)數(shù)據(jù)來(lái)修正更新系統(tǒng)狀態(tài)量的預(yù)測(cè)和其協(xié)方差的預(yù)測(cè)，從而獲得系統(tǒng)狀態(tài)量的估計(jì).所有的實(shí)際觀測(cè)值以迭代的方式都有助于系統(tǒng)狀態(tài)和原子干涉波形的估計(jì).船載絕對(duì)重力測(cè)量系統(tǒng)的觀測(cè)模型可描述為:

其中，y是原子布居數(shù)，A是條紋的直流偏置，C是條紋的對(duì)比度，keff是拉曼光有效波矢，可以表示為keff≈4π/λ，λ為激光波長(zhǎng)，g是絕對(duì)重力值，T是拉曼光脈沖的時(shí)間間隔，α是拉曼光的頻率啁啾率，vN是觀測(cè)噪聲.

考慮到α和φvib是比較準(zhǔn)確的輸入量，且連續(xù)兩次測(cè)量之間基本相互獨(dú)立.因此，未將它們作為系統(tǒng)隨時(shí)間變化的狀態(tài)參量.根據(jù)觀測(cè)模型，船載絕對(duì)重力測(cè)量系統(tǒng)可以用一個(gè)三參量的狀態(tài)向量X進(jìn)行描述，其形式如下:

狀態(tài)向量隨時(shí)間的演化關(guān)系可以用如下方程描述:

其中，Xk表示k時(shí)刻的系統(tǒng)狀態(tài)，Xk-1表示k-1 時(shí)刻的系統(tǒng)狀態(tài)，F(xiàn)是狀態(tài)向量X的演化矩陣，w是狀態(tài)演化過(guò)程中狀態(tài)向量X的過(guò)程噪聲.演化矩陣F可以通過(guò)狀態(tài)變量隨時(shí)間的演化關(guān)系獲得.由于這些隨機(jī)變量基本獨(dú)立，狀態(tài)向量X的過(guò)程噪聲w的協(xié)方差矩陣Q可以用一個(gè)對(duì)角矩陣表示.F和Q可表示為:

其中，σA，σC和σg分別是通過(guò)統(tǒng)計(jì)方法獲得的條紋直流偏置、對(duì)比度和重力值的標(biāo)準(zhǔn)偏差.

船載冷原子干涉式重力測(cè)量系統(tǒng)的重力值估計(jì)過(guò)程可以分為兩個(gè)步驟:狀態(tài)的預(yù)測(cè)和狀態(tài)的更新.首先，通過(guò)前一次測(cè)量的狀態(tài)估計(jì)對(duì)當(dāng)前時(shí)刻的狀態(tài)進(jìn)行預(yù)測(cè)，相應(yīng)的表達(dá)式如下:

其中，協(xié)方差矩陣P表示狀態(tài)預(yù)測(cè)的不確定性.是前一次測(cè)量的協(xié)方差矩陣的估計(jì)，是當(dāng)前測(cè)量協(xié)方差矩陣的預(yù)測(cè)估計(jì).其次，利用當(dāng)前的實(shí)際測(cè)量值對(duì)當(dāng)前時(shí)刻狀態(tài)和協(xié)方差的預(yù)測(cè)進(jìn)行更新，從而得到當(dāng)前時(shí)刻系統(tǒng)狀態(tài)的估計(jì)和協(xié)方差的估計(jì)，它們將作為下一時(shí)刻的預(yù)測(cè)初始量，其形式如下:

其中，Kk是卡爾曼濾波增益，它控制系統(tǒng)狀態(tài)量的更新，y() 是對(duì)當(dāng)前時(shí)刻觀測(cè)值的估計(jì)，yk是當(dāng)前時(shí)刻的實(shí)際觀測(cè)值，I是單位對(duì)角矩陣，H是觀測(cè)矩陣(狀態(tài)向量關(guān)于觀測(cè)方程的雅克比矩陣).Hk和Kk的形式如下:

其中，R是原子干涉條紋測(cè)量噪聲的方差，卡爾曼濾波增益Kk是估計(jì)值的不確定性與測(cè)量噪聲之間妥協(xié)的結(jié)果，EKF 算法中狀態(tài)量的噪聲參數(shù)可以通過(guò)先前的或?qū)嶋H的測(cè)量數(shù)據(jù)統(tǒng)計(jì)獲得.

為了消除與拉曼光波矢無(wú)關(guān)的相移的影響，船載絕對(duì)重力測(cè)量系統(tǒng)采用拉曼光正負(fù)掃頻的方式進(jìn)行動(dòng)態(tài)測(cè)量.數(shù)據(jù)處理時(shí)按時(shí)間測(cè)量順序使用(5)—(8)式分別對(duì)正負(fù)掃描的原子干涉條紋數(shù)據(jù)進(jìn)行迭代遞推處理，通過(guò)EKF 算法的狀態(tài)預(yù)測(cè)和更新，估計(jì)出每次下落測(cè)量時(shí)系統(tǒng)的狀態(tài)量X，其中g(shù)是估計(jì)的絕對(duì)重力值.將正負(fù)掃頻的重力測(cè)量值進(jìn)行交替平均處理，從而抵消與拉曼光波矢無(wú)關(guān)的相移所引入的重力值偏差.最后，進(jìn)行系統(tǒng)誤差的修正，獲得最終的絕對(duì)重力數(shù)據(jù).

3 船載絕對(duì)重力動(dòng)態(tài)測(cè)量實(shí)驗(yàn)系統(tǒng)

船載冷原子干涉式絕對(duì)重力動(dòng)態(tài)測(cè)量系統(tǒng)的示意如圖2 所示.原子重力儀、慣性穩(wěn)定平臺(tái)等儀器設(shè)備集成在由集裝箱改造而成的可移動(dòng)原子重力儀實(shí)驗(yàn)室內(nèi)，可移動(dòng)實(shí)驗(yàn)室安裝在試驗(yàn)船后部的甲板上，并與船只剛性連接，以應(yīng)對(duì)惡劣的海況條件，也可以保障儀器設(shè)備在良好的環(huán)境條件下正常工作運(yùn)轉(zhuǎn).

圖2 船載絕對(duì)重力動(dòng)態(tài)測(cè)量系統(tǒng)示意圖Fig.2.Schematic diagram of marine dynamic absolute gravity measurement system.

可移動(dòng)實(shí)驗(yàn)室包含的儀器設(shè)備有原子重力儀、慣性穩(wěn)定平臺(tái)、不間斷電源(uninterruptible power supply，UPS)供電系統(tǒng)、空調(diào)溫控裝置、主控電腦、差分定位系統(tǒng)(differential global position system，DGPS)、集裝箱等，詳細(xì)的布局及內(nèi)部照片如圖3所示.原子重力儀為本研究小組自行研發(fā)，成熟度、小型化程度和精度較高[5，12，14，27，28].慣性穩(wěn)定平臺(tái)的主要作用是實(shí)時(shí)維持重力傳感器的姿態(tài)，使其一直處于豎直方向，確保拉曼光方向與重力方向?qū)崟r(shí)重合.UPS 供電系統(tǒng)為可移動(dòng)實(shí)驗(yàn)室內(nèi)的所有儀器供電.差分定位系統(tǒng)安裝在集裝箱頂部，其定位精度可達(dá)cm 級(jí)別，可以記錄空間位置信息，包括經(jīng)緯度、海拔高度、船速和航向等.空調(diào)溫控裝置維持控制儀器的外部環(huán)境溫度和濕度，以確保儀器的正常工作.

圖3 可移動(dòng)的原子重力儀實(shí)驗(yàn)室 (a)實(shí)驗(yàn)室內(nèi)部?jī)x器設(shè)備布局示意圖;(b)可移動(dòng)實(shí)驗(yàn)室實(shí)物圖Fig.3.Transportable laboratory for atomic gravimeter:(a) Schematic diagram of the internal layout of instruments and equipments in the laboratory;(b) photo of the transportable laboratory.

4 船載絕對(duì)重力動(dòng)態(tài)測(cè)量實(shí)驗(yàn)結(jié)果

4.1 測(cè)量環(huán)境及流程

本文的絕對(duì)重力動(dòng)態(tài)測(cè)量實(shí)驗(yàn)場(chǎng)地位于南海某海域，測(cè)試時(shí)間在2020 年9 月末.在測(cè)量船航行前，所有的儀器設(shè)備都被調(diào)至最佳的運(yùn)行狀態(tài).此外，當(dāng)船停泊在碼頭時(shí)，通過(guò)測(cè)量不同T下的原子干涉條紋，初步評(píng)估了起始點(diǎn)的絕對(duì)重力值，并確定了動(dòng)態(tài)重力測(cè)量時(shí)合適的T和拉曼光頻率啁啾率的掃描范圍.此次動(dòng)態(tài)測(cè)量實(shí)驗(yàn)中，開(kāi)展了不同船速下的絕對(duì)重力測(cè)量以進(jìn)行比較.拉曼光脈沖的時(shí)間間隔T是4 ms，而原子團(tuán)下落的頻率是2 Hz.圖4(a)和(b)分別是動(dòng)態(tài)測(cè)量的航行路線和航行速度.其中，路線總長(zhǎng)度約85 km，船的最大速度約24.2 km/h，紅色和綠色航線上船航行的平均速度約22.6 km/h，相應(yīng)的路線長(zhǎng)度分別為42.9 km 和28.5 km;藍(lán)色航線所對(duì)應(yīng)的平均航行速度約14.3 km/h，相應(yīng)的路線長(zhǎng)度為12.8 km;黃色線路的航速基本小于2.1 km/h，這段航線包含了試驗(yàn)船在海上系泊和作業(yè)的狀態(tài)，船的活動(dòng)范圍未超過(guò)0.5 km.由于本次動(dòng)態(tài)試驗(yàn)主要以搭載的形式開(kāi)展，測(cè)量船以其他的科研任務(wù)為主，有固定的路線規(guī)劃，所以并未進(jìn)行重復(fù)測(cè)線和交叉測(cè)線等實(shí)驗(yàn).

圖4 絕對(duì)重力動(dòng)態(tài)測(cè)量的航線與船速 (a)航行路線;(b)航行速度Fig.4.Route and ship speed of absolute gravity dynamic measurement:(a) Sailing route of the ship;(b) speed of the ship.

海上測(cè)試期間，可移動(dòng)實(shí)驗(yàn)室外部的溫度和濕度分別在29 ℃和85%左右，相應(yīng)的變化幅度分別超過(guò)了6 ℃和30%.經(jīng)過(guò)可移動(dòng)實(shí)驗(yàn)室內(nèi)部的溫濕度控制后，溫度和濕度分別維持在20 ℃和55%左右.與靜態(tài)的環(huán)境相比，海上動(dòng)態(tài)測(cè)量時(shí)的振動(dòng)噪聲增大了幾個(gè)數(shù)量級(jí)，不同航速下振動(dòng)的功率譜密度(PSD)如圖5(a)所示，由于受海浪影響，振動(dòng)功率譜密度在0.23 Hz 附近有個(gè)比較大的特征峰.平均航速在22.6 km/h 的條件下，瞬時(shí)最大振動(dòng)加速度超過(guò)了0.158 m/s2.因此動(dòng)態(tài)測(cè)量時(shí)面臨的實(shí)驗(yàn)環(huán)境較為惡劣.

圖5 船載動(dòng)態(tài)測(cè)量的振動(dòng)環(huán)境與原子干涉條紋 (a)豎直方向振動(dòng)加速度的功率譜密度(PSD);(b)原子干涉條紋Fig.5.The vibration environment and atomic interference fringes for shipborne dynamic measurements:(a) Power spectral density(PSD) of the vertical vibration acceleration;(b) atomic interference fringes.

結(jié)合原子干涉儀的傳遞函數(shù)和振動(dòng)的噪聲功率譜密度[29]，可以計(jì)算動(dòng)態(tài)條件下振動(dòng)噪聲對(duì)原子干涉儀相位噪聲的影響.船速22.6 km/h，14.3 km/h和小于2.1 km/h 的條件下振動(dòng)對(duì)單次測(cè)量的相位噪聲貢獻(xiàn)分別可達(dá)102.1 rad，61.7 rad 和25.8 rad，振動(dòng)噪聲的貢獻(xiàn)遠(yuǎn)大于重力值隨時(shí)間和空間變化所引起的相位改變.圖5(b)是原子干涉條紋，其中，黑色數(shù)據(jù)點(diǎn)是實(shí)際測(cè)量的原子干涉條紋數(shù)據(jù)，難以分辨出有效的原子干涉波形.將振動(dòng)加速度所產(chǎn)生的干擾相位φvib進(jìn)行修正后[16]，還原出了較為有效的原子干涉條紋(紅色菱形數(shù)據(jù)點(diǎn)).但由于加速計(jì)性能的限制，并不能完全獲得拉曼反射鏡真實(shí)的振動(dòng)情況，還存在一定的殘余相位噪聲干擾，測(cè)量數(shù)據(jù)的噪聲水平依舊較大.因此，本文使用EKF 算法對(duì)測(cè)量數(shù)據(jù)進(jìn)行濾波迭代處理，進(jìn)一步從噪聲中分離出有效的絕對(duì)重力值.

4.2 基于擴(kuò)展卡爾曼濾波算法的絕對(duì)重力測(cè)量數(shù)據(jù)處理結(jié)果

基于第2 小節(jié)原理部分介紹的數(shù)據(jù)處理方法，獲得了由EKF 算法估計(jì)的絕對(duì)重力值數(shù)據(jù)gEst(紅色曲線).在原子干涉條紋還原的基礎(chǔ)上，使用最小二乘法的余弦擬合獲得了的絕對(duì)重力數(shù)據(jù)gCorr(黑色曲線).通過(guò)固體潮汐模型對(duì)重力數(shù)據(jù)進(jìn)行了修正，相應(yīng)結(jié)果如圖6 所示，gEst和gCorr對(duì)應(yīng)的采樣率分別為2 Hz 和1/59 Hz.利用這些數(shù)據(jù)分析了船載絕對(duì)重力動(dòng)態(tài)測(cè)量的靈敏度和準(zhǔn)確度，下面分別進(jìn)行介紹.

圖6 絕對(duì)重力值原始數(shù)據(jù)Fig.6.Raw data for absolute gravity values.

4.2.1 船載絕對(duì)重力動(dòng)態(tài)測(cè)量靈敏度分析

在船速小于2.1 km/h 的條件下，原子干涉條紋數(shù)據(jù)的噪聲相對(duì)較小且比較穩(wěn)定，利用相應(yīng)的重力數(shù)據(jù)可以計(jì)算Allan 偏差，以評(píng)估動(dòng)態(tài)條件下重力測(cè)量的靈敏度和噪聲水平，相應(yīng)的Allan 偏差分析結(jié)果如圖7 所示.其中，藍(lán)色圓點(diǎn)是gCorr的Allan偏數(shù)據(jù)，相應(yīng)的藍(lán)色曲線是通過(guò)τ-1/2擬合所得的結(jié)果，對(duì)應(yīng)的靈敏度為300.2 mGal/Hz1/2;紅色方塊是通過(guò)gEst分析的Allan 偏差數(shù)據(jù)，紅色曲線是其相應(yīng)的擬合結(jié)果，對(duì)應(yīng)的靈敏度約為136.8 mGal/Hz1/2.因此，用EKF 算法估計(jì)的重力值具備更低的噪聲水平，靈敏度比最小二乘法的結(jié)果提高了約53.8%，EKF 算法對(duì)動(dòng)態(tài)測(cè)量下的測(cè)量噪聲有較好的抑制效果.Allan 偏差數(shù)據(jù)可以反映重力測(cè)量值的高斯白噪聲水平，EKF 算法在給定的噪聲參數(shù)條件下，相當(dāng)于具備100 s 低通濾波的效果，同時(shí)，由于惡劣動(dòng)態(tài)環(huán)境下系統(tǒng)穩(wěn)定性的原因，重力測(cè)量值存在有色噪聲，這部分噪聲通過(guò)算法無(wú)法被很好處理，因此100 s 積分時(shí)間以下的數(shù)據(jù)存在由低到高的趨勢(shì)，且低于理想的擬合曲線.

圖7 通過(guò)不同算法獲得的重力值的Allan 偏差數(shù)據(jù)Fig.7.Allan deviation data of gravity values obtained by different algorithms.

在EKF 算法中，需要使用系統(tǒng)狀態(tài)量的方差和原子布居數(shù)測(cè)量噪聲的方差作為其參數(shù)，這些通過(guò)先驗(yàn)的實(shí)驗(yàn)數(shù)據(jù)進(jìn)行統(tǒng)計(jì)獲得.其中，測(cè)量噪聲的方差包含了探測(cè)噪聲和振動(dòng)相位修正后殘余相位噪聲引入的原子布居數(shù)波動(dòng).探測(cè)噪聲引入的原子布居數(shù)噪聲的方差可以通過(guò)文獻(xiàn)[11]中給出的原子布居數(shù)的統(tǒng)計(jì)概率公式進(jìn)行分析和評(píng)估.拉曼光負(fù)掃頻時(shí)，由原子布居數(shù)統(tǒng)計(jì)出的偏置、對(duì)比度和探測(cè)噪聲引起的幅度噪聲分別為0.502±0.007，0.109±0.004 和0.036±0.007，正掃頻 時(shí)分別 為0.482±0.008，0.128±0.002 和0.035±0.007.除 了需要考慮探測(cè)噪聲對(duì)原子布居數(shù)噪聲的影響外，還需評(píng)估殘余振動(dòng)噪聲以及其他噪聲的影響.通過(guò)最小二乘法余弦擬合獲得的重力數(shù)據(jù)包含了所有噪聲的影響.當(dāng)剔除探測(cè)引起的原子布居噪聲的影響后，大致評(píng)估出剩余相位噪聲的標(biāo)準(zhǔn)偏差為0.744 rad.考慮分析出的探測(cè)和剩余相位噪聲的影響，通過(guò)理想條件下的模擬計(jì)算，可以得出其對(duì)原子布居數(shù)噪聲的總貢獻(xiàn)，相應(yīng)的標(biāo)準(zhǔn)偏差為0.06，方差是0.0036.將其作為EKF 算法中測(cè)量噪聲的方差參數(shù)，這相當(dāng)于減小了重力值估計(jì)的不確定性.

在單次測(cè)量中，空間位置和時(shí)間變化引入的重力值變化在1 mGal 以?xún)?nèi)，其不確定性相對(duì)比較小.但是，為了兼顧算法中重力值迭代的收斂速度，其相應(yīng)的方差使用了實(shí)驗(yàn)室條件下對(duì)應(yīng)T的重力測(cè)量值方差，通過(guò)若干觀測(cè)數(shù)據(jù)就可以收斂.在初始數(shù)據(jù)對(duì)整體重力值數(shù)據(jù)的影響較小時(shí)，可以進(jìn)一步優(yōu)化算法中重力值的方差參數(shù)，從而估計(jì)出不確定性更小的重力值數(shù)據(jù).通過(guò)測(cè)量噪聲和系統(tǒng)狀態(tài)量不確定性之間的權(quán)衡，達(dá)到抑制測(cè)量噪聲和提升靈敏度的目的.

4.2.2 船載絕對(duì)重力動(dòng)態(tài)測(cè)量的準(zhǔn)確度分析

將不同船速下的測(cè)量數(shù)據(jù)處理后，獲得了不同算法計(jì)算的絕對(duì)重力數(shù)據(jù)gCorr(基于修正振動(dòng)相位干擾后的原子干涉條紋數(shù)據(jù)，通過(guò)余弦擬合所得的絕對(duì)重力值，采樣周期59 s)、gEKF(基于EKF算法濾波估計(jì)的原子干涉條紋數(shù)據(jù)，通過(guò)余弦擬合所得的絕對(duì)重力值，采樣周期59 s)和gEst(基于EKF 算法狀態(tài)估計(jì)所得的絕對(duì)重力值，采樣周期0.5 s)，這些重力值的對(duì)比結(jié)果如圖8 所示，單個(gè)數(shù)據(jù)點(diǎn)經(jīng)過(guò)5 min 時(shí)間的平均，并減去了同一絕對(duì)重力參考值g0.圖8(a)是不同計(jì)算方法獲得的未修正系統(tǒng)誤差的絕對(duì)重力值對(duì)比結(jié)果;圖8(b)是修正系統(tǒng)誤差后的絕對(duì)重力值與重力模型XGM2019計(jì)算結(jié)果的對(duì)比;圖8(c)是絕對(duì)重力值gEst與重力模型的殘差.由圖8(a)知，3 種處理方式獲得的絕對(duì)重力值的變化趨勢(shì)基本一致，gCorr的噪聲水平最大，而gEKF和gEst的離散程度較小，并且兩者基本吻合.與最小二乘法的余弦擬合相比，使用EKF 算法可以有效的抑制部分噪聲，并且還可以通過(guò)觀測(cè)數(shù)據(jù)有效地追蹤估計(jì)動(dòng)態(tài)條件下的重力值;同時(shí)，迭代的方式使采樣率更高，從而可以較好地反映重力值隨時(shí)間的變化.

圖8 動(dòng)態(tài)絕對(duì)重力數(shù)據(jù)處理結(jié)果 (a)使用不同算法獲得的絕對(duì)重力值;(b)修正系統(tǒng)誤差后的絕對(duì)重力值;(c)絕對(duì)重力值gEst 與由重力模型計(jì)算的結(jié)果的殘差數(shù)據(jù)Fig.8.Comparison of the absolute gravity data:(a) Absolute gravity values obtained by the different algorithms;(b)absolute gravity values obtained after correcting the systematic errors;(c) the residual data between the absolute gravity values gEst and the calculated results based on the gravity model.

圖8(b)是絕對(duì)重力值gEKF和gEst在修正各項(xiàng)系統(tǒng)誤差后與絕對(duì)重力模型XGM2019 計(jì)算的結(jié)果進(jìn)行對(duì)比的情況.重力模型XGM2019 是一個(gè)組合的全球重力場(chǎng)模型，由度和階(d/o)高達(dá)5399的球面諧波表示，對(duì)應(yīng)于2'(約4 km)的空間分辨率，它使用衛(wèi)星數(shù)據(jù)、地面重力網(wǎng)等多種數(shù)據(jù)源作為參考，與之前的模型(如XGM2016，EIGEN6c4或EGM2008)相比部分性能更優(yōu)[30].絕對(duì)重力數(shù)據(jù)與由重力模型計(jì)算的絕對(duì)重力值基本一致，且比較吻合.

一般靜態(tài)條件下，原子重力儀的系統(tǒng)誤差主要包括雙光子光移、科里奧利力效應(yīng)、拉曼光波前畸變等[5，12，14，27，31]，除了雙光子光移外，系統(tǒng)誤差的修正量為(41.2±11.4) μGal.由于動(dòng)態(tài)條件下T較小，雙光子效應(yīng)根據(jù)理論計(jì)算[32]，評(píng)估的修正量約為(—2.738±0.262) mGal.動(dòng)態(tài)航行條件下，地球自轉(zhuǎn)和船體運(yùn)動(dòng)會(huì)產(chǎn)生額外的垂向加速度干擾，即E?tv?s 效應(yīng)[17，18，33]，該效應(yīng)的項(xiàng)修量在—45.661—39.516 mGal 之間，相應(yīng)的不確定度為0.343 mGal.經(jīng)評(píng)估，船載原子絕對(duì)重力測(cè)量系統(tǒng)的總系統(tǒng)誤差修正量基本介于—48.440—36.737 mGal之間，相應(yīng)的不確定度為0.434 mGal (k=1).

將絕對(duì)重力值gEst與重力模型計(jì)算的絕對(duì)重力值相減后得到了殘差數(shù)據(jù)，如圖8(c)所示，全部重力殘差數(shù)據(jù)的均值為(0.8±10.8) mGal (k=1).通過(guò)紅色虛線標(biāo)記的區(qū)間1 和4 所對(duì)應(yīng)的平均船速約為22.6 km/h，相應(yīng)的殘差均值為(2.6±13.6)mGal (k=1);區(qū)間2 的平均船速是14.3 km/h，相應(yīng)的殘差均值為(3.4±15.3) mGal (k=1);區(qū)間3 的船速小于2.1 km/h，相應(yīng)的殘差均值為(—0.3±8.1) mGal (k=1).重力模型XGM2019 在海面、平原等地形平坦區(qū)具有相對(duì)較高的精度，一般在幾毫伽以?xún)?nèi).與重力模型的比較粗略反映了測(cè)量精度，更精確的評(píng)估需要后續(xù)更多的動(dòng)態(tài)測(cè)量實(shí)驗(yàn)來(lái)驗(yàn)證，如與傳統(tǒng)海洋重力儀的比對(duì)、交叉測(cè)點(diǎn)和重復(fù)測(cè)線等實(shí)驗(yàn)，后續(xù)實(shí)驗(yàn)仍在開(kāi)展之中.相關(guān)的噪聲也需要進(jìn)一步分析，從而進(jìn)行改進(jìn)以提升系統(tǒng)的整體測(cè)量性能.

目前，限制船載絕對(duì)重力測(cè)量系統(tǒng)性能的主要因素有振動(dòng)噪聲、旋轉(zhuǎn)噪聲和探測(cè)噪聲等.由于加速度計(jì)的測(cè)量精度、靈敏度有限，尚存在一定的殘余噪聲，因此需要使用更高精度、更高性能的加速度計(jì)采集拉曼光反射鏡的振動(dòng)加速度信號(hào).為了維持重力傳感器的水平狀態(tài)，慣性穩(wěn)定平臺(tái)的姿態(tài)需要實(shí)時(shí)調(diào)節(jié)，實(shí)時(shí)的隨機(jī)旋轉(zhuǎn)也會(huì)引入測(cè)量噪聲，這需要通過(guò)合適的評(píng)估方法來(lái)處理這一部分噪聲干擾.由于測(cè)量環(huán)境惡劣，原子數(shù)探測(cè)和計(jì)算的準(zhǔn)確性也會(huì)引入一定的隨機(jī)噪聲，需要通過(guò)軟硬件結(jié)合的方式進(jìn)行改進(jìn).算法上，振動(dòng)噪聲和探測(cè)噪聲以外的其他殘余噪聲的影響沒(méi)有較好的關(guān)系式或者統(tǒng)計(jì)上的方法進(jìn)行分析和描述，無(wú)法很好地將它們納入到狀態(tài)模型中，只能進(jìn)行整體上的考量，這也是限制靈敏度提升的主要因素，同時(shí)也是改進(jìn)方向之一.上述內(nèi)容也是接下來(lái)提升船載絕對(duì)重力測(cè)量系統(tǒng)性能的后續(xù)研究方向.

5 結(jié)論

本文基于自主研制的船載冷原子干涉式絕對(duì)重力測(cè)量系統(tǒng)，實(shí)現(xiàn)了海上不同航速條件下的動(dòng)態(tài)絕對(duì)重力測(cè)量，研究了基于EKF 的動(dòng)態(tài)數(shù)據(jù)處理方法.為了抑制原子干涉條紋數(shù)據(jù)的殘余噪聲，本文從算法角度出發(fā)引入了適用于該系統(tǒng)的EKF 算法進(jìn)行數(shù)據(jù)的處理.與最小二乘余弦擬合的方法相比，EKF 算法可以將測(cè)量靈敏度從300.2 mGal/Hz1/2提升到136.8 mGal/Hz1/2(T=4 ms).因此，EKF 算法在使用合理評(píng)估的狀態(tài)量噪聲參數(shù)的條件下對(duì)重力測(cè)量的噪聲有較好的抑制效果.

EKF 算法通過(guò)觀測(cè)數(shù)據(jù)可以有效地追蹤估計(jì)動(dòng)態(tài)條件下的重力值.同時(shí)，以迭代的方式使采樣率更高，實(shí)時(shí)性更強(qiáng).通過(guò)EKF 算法估計(jì)的絕對(duì)重力值的噪聲水平較小，而且與最小二乘法的余弦擬合結(jié)果的變化趨勢(shì)基本一致，證實(shí)了EKF 算法對(duì)重力值的估計(jì)是可靠和準(zhǔn)確的.將修正各項(xiàng)系統(tǒng)效應(yīng)后的絕對(duì)重力數(shù)據(jù)與由重力場(chǎng)模型XGM2019計(jì)算的結(jié)果進(jìn)行了比較，獲得了航速22.6 km/h，14.3 km/h 和小于2.1 km/h 條件下殘差的均值，分別為(2.6±13.6) mGal，(3.4±15.3) mGal 和(—0.3±8.1) mGal，整體上所以殘差數(shù)據(jù)的均值為(0.8±10.8) mGal (k=1).本文將EKF 作為船載原子絕對(duì)重力測(cè)量系統(tǒng)的估計(jì)器，為重力值的估計(jì)和重力測(cè)量噪聲的抑制提供了一種有效的處理方法和新思路.