李 達，張海洋，李城鎖，向民志，高 巍，李 中,5

（1. 哈爾濱工程大學(xué) 智能科學(xué)與工程學(xué)院，哈爾濱 150001；2. 天津航海儀器研究所，天津 300131；3. 近地面探測技術(shù)重點實驗室，無錫 214035；4. 信息工程大學(xué) 地理空間信息學(xué)院，鄭州 450001；5. 中船集團航海保障技術(shù)實驗室，天津 300131）

重力梯度是重力矢量的空間變化率，在礦產(chǎn)資源勘探、地球科學(xué)研究和慣性導(dǎo)航等領(lǐng)域具有重要意義[1,2]。重力梯度儀是用于測量重力梯度的精密設(shè)備，Bell Aerospace公司提出的基于旋轉(zhuǎn)加速度計測量原理的重力梯度儀是迄今唯一實用的近地表動態(tài)重力梯度儀[3]。它通過對稱反向安裝的加速度計抵御載體水平線運動對重力梯度測量的影響，然而地表附近真實重力梯度異常微弱，通常在十幾到幾百E量級（1 E=10-9s-2），而10 E的重力梯度在相距100 mm的兩點只能產(chǎn)生約為10-10g的引力加速度差異，這對于核心敏感元件加速度計與儀器動態(tài)適應(yīng)性構(gòu)成了嚴峻挑戰(zhàn)[4]。旋轉(zhuǎn)加速度計式重力梯度儀通過機械旋轉(zhuǎn)的方式將重力梯度信號調(diào)制到敏感器輸出信號旋轉(zhuǎn)頻率的二倍頻分量上，有效降低了對加速度計零位漂移與測量分辨率的要求。而在動態(tài)重力梯度測量中，對稱反向安裝的加速度計間微小的標度因數(shù)差異會使得載體水平線加速度進入重力梯度敏感器輸出信號，形成測量誤差[5]。為實現(xiàn)高精度重力梯度動態(tài)測量，要求對稱反向安裝的兩組加速度計間標度因數(shù)不一致性優(yōu)于0.1 ppm以內(nèi)[6]，這對加速度計間標度因數(shù)不一致程度提出了嚴苛的要求。

為此，諸多研究人員提出了加速度計標度因數(shù)調(diào)整方案[7,8]，分別采用不同的調(diào)整方案取得了較好的控制效果，但絕大多數(shù)控制方案只是針對靜態(tài)下的加速度計標度因數(shù)調(diào)整，而在動態(tài)下如何應(yīng)對載體運動干擾的問題缺乏深入研究。尤其在航空測量中，受空中橫風(fēng)及大氣湍流影響，搭載飛機經(jīng)常出現(xiàn)線加速度短期波動，影響調(diào)整回路工作。針對上述問題，本文提出基于Sage-Husa自適應(yīng)濾波的加速度計標度因數(shù)調(diào)整方法，實時調(diào)整濾波參數(shù)以適應(yīng)不同測量工況下的載體運動干擾，提高加速度計標度因數(shù)一致性調(diào)整精度。

1 誤差機理與標度因數(shù)調(diào)整回路

1.1 測量原理

如圖1所示，旋轉(zhuǎn)加速度計式重力梯度儀主體儀器由重力梯度敏感器和慣性穩(wěn)定平臺兩個關(guān)鍵組件構(gòu)成。重力梯度敏感器用于完成重力梯度張量水平分量的測量；慣性穩(wěn)定平臺則用于承載重力梯度敏感器，為動態(tài)重力梯度測量提供穩(wěn)定的動力學(xué)環(huán)境，抑制載體角運動對測量的影響，同時將重力梯度敏感器穩(wěn)定在地理坐標系下。如圖2所示，重力梯度敏感器基于加速度計位置差分的測量原理，通過機械旋轉(zhuǎn)的方式將旋轉(zhuǎn)中心處的重力梯度張量水平分量調(diào)制到系統(tǒng)旋轉(zhuǎn)頻率的二倍頻處，重力梯度敏感器輸出信號與重力梯度張量水平分量之間的關(guān)系可表示為[9]：

圖1 旋轉(zhuǎn)加速度計式重力梯度儀主體儀器構(gòu)成示意圖Fig.1 Schematic diagram of measurement principle of gravity gradient sensor

圖2 重力梯度敏感器測量原理示意圖Fig.2 Schematic diagram of measurement principle of gravity gradient sensor

式中，ai（i= 1,2,3,4）是第i只加速度計敏感軸方向的比力，l是加速度計檢測質(zhì)心到旋轉(zhuǎn)中心的距離，Γuv、Γxy是對應(yīng)方向的重力梯度張量分量（其中單位是E（1 E=10-9s-2），ω是旋轉(zhuǎn)圓盤的旋轉(zhuǎn)角速率。進行動態(tài)測量時，慣性穩(wěn)定平臺采用三環(huán)固定指北半解析式控制方案，在隔離載體角運動的同時，將重力梯度敏感器穩(wěn)定在地理坐標系，為重力梯度測量提供測量方向基準。

1.2 加速度計標度因數(shù)調(diào)整回路

針對高精度重力梯度測量需要兩組加速度計標度因數(shù)不一致性達到0.1 ppm量級的要求，Metzger E H首次提出了加速度計標度因數(shù)一致性調(diào)整方案[10]，具體實施思路為：（1）在加速度計表體磁路中設(shè)置調(diào)整機構(gòu)，通過控制標度調(diào)整電流實現(xiàn)加速度計標度因數(shù)在線調(diào)整功能；（2）由技術(shù)設(shè)計和工藝控制保證配對的加速度計標度因數(shù)初始差異不大于1‰；（3）傾斜敏感器使重力分量在敏感器X方向形成常值分量，激勵兩組加速度計標度因數(shù)不一致信息調(diào)制在旋轉(zhuǎn)頻率的正余弦分量上，再通過同步解調(diào)提??；（4）對標度因數(shù)不一致信息進行PID校正，控制加速度計表體調(diào)整電流，實現(xiàn)加速度計標度因數(shù)實時調(diào)整功能。

靜態(tài)條件下，加速度計標度因數(shù)的變化緩慢，因此容易通過調(diào)整PID校正環(huán)節(jié)的參數(shù)實現(xiàn)3#、4#加速度計標度因數(shù)跟蹤1#、2#加速度計。而在動態(tài)條件下，載體水平加速度是時變的，通過解調(diào)提取的兩路加速度計標度因數(shù)不一致信息存在動態(tài)誤差，導(dǎo)致PID校正環(huán)節(jié)輸出的調(diào)整電流存在波動，形成調(diào)整誤差。在重力梯度動態(tài)測量中，標度因數(shù)不一致誤差引起水平加速度進入重力梯度測量通道，已逐漸成為一項重要誤差來源。

2 加速度計標度因數(shù)動態(tài)實時調(diào)整方法

2.1 Kalman濾波方程

為提高動態(tài)下兩路加速度計標度因數(shù)不一致信息提取精度，在不一致信息提取中引入Kalman濾波環(huán)節(jié)，對兩路信息實現(xiàn)最優(yōu)估計。系統(tǒng)的狀態(tài)轉(zhuǎn)移方程離散形式為：

式中，X=[ ΔK13ΔK24]T，ΔK13和ΔK24為兩路加速度計標度因數(shù)不一致數(shù)值，βi（i= 1,2,3,4）是事先標定的第i只加速度計標度因數(shù)溫度系數(shù)，Δti（i= 1,2,3,4）是第i只加速度計的采集溫度變化量，W是2維系統(tǒng)噪聲向量。

以重力梯度儀一倍頻解調(diào)值作為觀測量，系統(tǒng)的量測方程離散形式為：

式中Z=[C1C2]T，C1和C2是經(jīng)解調(diào)提取的重力梯度敏感器輸出信號旋轉(zhuǎn)頻率正余弦分量，V是2維量測噪聲向量，量測矩陣H為：

在Kalman濾波參數(shù)中，系統(tǒng)噪聲方差陣Q的數(shù)值可通過重力梯度儀靜態(tài)數(shù)據(jù)計算得到，但量測噪聲方差陣R則與載體動態(tài)息息相關(guān)，很難通過事先預(yù)設(shè)的數(shù)值滿足同一種載體不同工況下的濾波要求。而R陣的設(shè)置誤差會導(dǎo)致Kalman濾波的精度降低，嚴重時還可能會導(dǎo)致濾波發(fā)散，在實際工程應(yīng)用中必須予以考慮。

2.2 Sage-Husa自適應(yīng)濾波

為解決上述問題，引入Sage-Husa自適應(yīng)濾波的方法估計R陣的數(shù)值[11,12]。定義量測預(yù)測誤差為：

由此得到R陣的表達式為：

為提高系統(tǒng)實時運算效率，將式(7)寫成遞推估計的形式，即：

當(dāng)k→∞時有即長時間濾波后自適應(yīng)能力將逐漸減弱，直至幾乎失去自適應(yīng)能力，為始終保持R陣的自適應(yīng)能力，將等加權(quán)平均改為指數(shù)漸消記憶加權(quán)平均，即：

式中初值λ0=1，c是漸消因子，要求0＜c＜1。當(dāng)k→∞時有λk→1-c，使濾波器始終保持R陣的自適應(yīng)能力。此外，為保證R陣的正定性，采用序貫濾波的方法對R陣的對角線元素進行限制，簡記則：

圖3 Sage-Husa自適應(yīng)濾波算法流程Fig.3 Sage-Husa adaptive filtering algorithm flow

2.3 漸消因子數(shù)值確定

在Sage-Husa自適應(yīng)濾波算法中，漸消因子c的參數(shù)選擇尤為重要。漸消因子c的取值越小，對最新量測噪聲變化的自適應(yīng)能力越強，但若漸消因子c的取值過小，噪聲方差陣的估計結(jié)果會變化劇烈，引起Kalman濾波的精度降低。為此，漸消因子c的取值范圍通常在0.9和0.999之間。本文采用半實物仿真的方式選擇漸消因子c的最佳參數(shù)。利用靜態(tài)下得到兩組加速度計標度因數(shù)不一致的變化過程，以及重力梯度儀航空試驗中的實測載體線加速度，結(jié)合重力梯度測量誤差方程建立重力梯度敏感器模擬輸出信號asum，仿真計算頻率為1 Hz。使用不同參數(shù)的濾波方法得到的兩組加速度計標度因數(shù)不一致信號如圖4所示。為評價不同參數(shù)的優(yōu)劣，定義提取誤差ε為：

式中，ε是本次仿真定義的提取誤差，1S和S2是兩路加速度計標度因數(shù)差異真值，和是兩路加速度計標度因數(shù)差異估計值，n是數(shù)據(jù)總數(shù)。

對同一組仿真數(shù)據(jù)使用常規(guī)解調(diào)、Kalman濾波以及不同漸消因子的Sage-Husa自適應(yīng)濾波的方法提取兩路加速度計標度因數(shù)不一致信號，與真值的對比如圖4所示，利用式(12)的計算方法統(tǒng)計不同信號處理方法的提取誤差，相關(guān)結(jié)果如表1所示。

結(jié)合圖4和表1可知，當(dāng)漸消因子c取0.95時，可實現(xiàn)兩路提取誤差ε最小且小于0.1 ppm，因此在本研究中漸消因子c的值取0.95。

圖4 不同提取方法的加速度計標度因數(shù)不一致信號提取結(jié)果Fig.4 Results of accelerometer scale factor inconsistent signal extraction with different extraction methods

表1 不同提取方法的加速度計標度因數(shù)不一致信號提取誤差統(tǒng)計Tab.1 Error statistics of accelerometer scale factor inconsistent signal extraction with different extraction methods

3 試驗驗證

3.1 加速度計內(nèi)部調(diào)整電流

動態(tài)條件下敏感器測量信號受到載體運動干擾，無法從其旋轉(zhuǎn)頻率正余弦分量準確獲得兩路加速度計標度因數(shù)不一致性信號，只能通過加速度計標度調(diào)整電流波動反推動態(tài)條件下加速度計標度因數(shù)調(diào)整精度。

項目研究團隊分別于2018年和2020年開展兩次基于旋轉(zhuǎn)加速度計原理的航空重力梯度測量試驗，首次航空試驗中加速度計標度因數(shù)控制回路采用常規(guī)解調(diào)的控制方案，第二次航空試驗采用基于Sage-Husa自適應(yīng)濾波的控制方案，兩次試驗中加速度計內(nèi)部調(diào)整電流波動情況分別如圖5和圖6所示。

圖5 首次航空試驗加速度計內(nèi)部調(diào)整電流波動Fig.5 Internal adjustment current fluctuation of accelerometer in the first aviation test

圖6 第二次航空試驗加速度計內(nèi)部調(diào)整電流波動Fig.6 Internal adjustment current fluctuation of accelerometer in the second aviation test

對比圖5和圖6，雖然加速度計標度調(diào)整電流無法體現(xiàn)加速度計標度因數(shù)不一致性調(diào)整回路精度，但對比標度調(diào)整電流在測量過程中的波動情況，第二次試驗中電流波動僅為第一次試驗的1/3，表明采用Sage-Husa自適應(yīng)濾波控制方案能夠抑制航空載體運動對加速度計標度因數(shù)調(diào)整的影響。

3.2 重力梯度測量結(jié)果及精度統(tǒng)計

兩次航空試驗中兩路重力梯度測量信號重復(fù)線測量典型結(jié)果如圖7和圖8所示，需要說明的是，由于兩次航空試驗測量地點不同，因此兩路重力梯度信號測量結(jié)果也不相同。

圖7 首次航空試驗兩路重力梯度信號處理結(jié)果Fig.7 Processing results of two gravity gradient signals in the first aviation test

圖8 第2次航空試驗兩路重力梯度信號處理結(jié)果Fig.8 Processing results of two gravity gradient signals in the second aviation test

鑒于目前國內(nèi)沒有高精度航空重力梯度儀進行外符合比對，只能采用內(nèi)符合的方式進行精度評價，內(nèi)符合精度計算公式為：

式中，εj是重力梯度j分量的內(nèi)符合精度，Γj1（i）和Γj2（i）是第1、2航次重力梯度j分量第i點的測量值，n是測量總點數(shù)。兩次航空試驗各條測量重復(fù)線內(nèi)符合精度統(tǒng)計結(jié)果如表2所示，結(jié)果表明調(diào)整回路采用Sage-Husa自適應(yīng)濾波控制方案時能夠大幅度提高航空重力梯度測量精度?？梢姡瑑?nèi)符合精度由150E@1 km提升至65E@1 km。

表2 兩次航空試驗內(nèi)符合精度統(tǒng)計Tab.2 Accuracy statistics in two aviation tests

4 結(jié) 論

本文針對旋轉(zhuǎn)加速度計式重力梯度儀在動基座測量過程中，加速度計標度因數(shù)不一致性調(diào)整回路易受到載體水平運動干擾且載體水平運動存在較大波動的問題，提出了一種Sage-Husa自適應(yīng)濾波的方法，實時適應(yīng)不同量級的載體運動干擾，對加速度計標度因數(shù)不一致信息進行最優(yōu)估計，并以此實時調(diào)整加速度計內(nèi)部電磁線圈電流來實現(xiàn)加速度計標度因數(shù)一致性調(diào)整。兩次航空試驗測量結(jié)果表明，較常規(guī)的解調(diào)方法，本文所提方法在相同工況下能夠?qū)⒓铀俣扔嫎硕日{(diào)整電流波動降低 2/3，測量內(nèi)符合精度由150E@1 km提高至65E@1km，證明了所提方法的有效性。