水平張量重力梯度儀垂向運動誤差實時補償方法

李 達，趙 明，范士鋒，李 中，李城鎖，趙 琳

（1.哈爾濱工程大學 智能科學與工程學院，哈爾濱 150001；2.天津航海儀器研究所，天津 300131；3.中國人民解放軍海軍裝備部，北京 100071；4.中國船舶航海保障技術實驗室，天津 300131）

重力梯度信號是重力矢量的空間微分信號，與傳統的重力異常相比，該信號擁有更高分辨率的地球重力場信息，在油氣資源勘探、地球內部構造研究以及慣性重力組合導航等領域具有重要意義[1-3]。目前存在靜電懸浮、冷原子干涉的超導等多種不同原理的重力梯度測量方案，但迄今僅有旋轉加速度計原理的重力梯度測量方案實現了近地表動態測量工程應用[4]。按照測量分量的不同，該原理的重力梯度儀分為水平張量式和全張量式兩種[5,6]。本文研究的重力梯度儀為水平張量式，能夠同時測量Γuv和Γxy兩種水平重力梯度張量分量。該原理的重力梯度儀將高分辨力加速度計敏感軸水平布放，用于抑制動態條件下載體垂向運動加速度產生的動態測量誤差。但由于加速度計安裝誤差的存在，載體垂向運動會以一定比例進入重力梯度測量通道，而真實的重力梯度信號極其微小。為此，一方面通過精密加工裝配盡可能減小安裝誤差[7]，另一方面，需要通過數據處理進一步抑制該項誤差對重力梯度測量的影響[8]。

在重力梯度垂向運動測量誤差數據處理方面，國外成熟商用重力梯度儀FALCONTM和Air-FTGTM的研究報告中，將垂向運動補償作為重力梯度數據處理的關鍵技術列出[9]，但詳細的數據處理方法并不對外公布。國內通常采用遞推最小二乘法估計加速度計敏感軸垂向角度綜合安裝誤差角，以此結合載體實測垂向加速度完成垂向運動誤差補償。該方法在短時間、海況平穩的情況下能夠較好地完成垂向運動誤差補償，但隨著測量時間的增加或者海況變化時，該方法對綜合安裝誤差角的估計精度變差，從而影響補償效果。本文根據船載實測重力梯度數據特征提出含有遺忘因子且具備自適應能力的垂向運動誤差實時補償方法，通過調整補償參數的方法能夠在不同測量工況下完成重力梯度垂向運動誤差的實時補償，為實現高精度動態重力梯度測量提供了一種數據實時處理方法。

1 垂向運動誤差機理分析

1.1 測量原理

重力梯度儀樣機由重力梯度敏感器和穩定平臺組成，并配套相應的減震器、環境控制裝置和顯控機箱等輔助測量設備，其中敏感器與慣性穩定平臺如圖1所示，敏感器用于測量水平分量重力梯度信號，穩定平臺則為其隔離載體角運動，并在動態環境下提供可靠的指向基準。重力梯度敏感器將四只高分辨力加速度計均勻布設在一個旋轉圓盤上，對稱安裝的加速度計敏感軸指向相反，用于減弱載體線運動對重力梯度測量的影響，正交安裝的加速度計則可以抑制旋轉誤差對重力梯度測量的影響，其原理圖如圖2 所示，測量原理公式[10]為：

圖1 重力梯度儀樣機組成Fig.1 Composition of gravity gradiometer

圖2 重力梯度敏感器原理圖Fig.2 Schematic diagram of gravity gradient sensor

式(1)中，ai(i=1,2,3,4)是第i只加速度計輸出的加速度測量值，l是加速度計旋轉半徑，Γuv、Γxy是重力梯度張量水平分量，其單位為E（1E=10-9s-2），ω是旋轉圓盤的旋轉角速度，一般為

1.2 誤差機理分析

重力梯度敏感器中的各加速度計敏感軸垂向角度安裝誤差如圖3 所示。

圖3 加速度計敏感軸垂向角度安裝誤差示意圖Fig.3 Schematic diagram of vertical angle installation error of accelerometer sensitive axis

該安裝誤差會使動態條件下的載體垂向運動干擾重力梯度敏感器測量信號，記asum=(a1+a3) -(a2+a4)，則有：

式(2)中δasum,az是由載體垂向運動引起的重力梯度敏感器測量誤差，az是載體垂向線加速度，βi(i=1,2,3,4)是第i只加速度計敏感軸的垂向角度安裝誤差。結合式(1)和式(2)可知，載體垂向運動中的2ω頻率分量，與綜合安裝誤差角(β1+β3) -(β2+β4)相耦合，形成重力梯度測量誤差。因此在重力梯度數據處理中，需要利用重力梯度儀慣性穩定平臺中垂向加速度計測量得到的載體垂向加速度，結合綜合安裝誤差角，進行載體垂向運動誤差補償。對式(2)進行全微分，同時記βsum=(β1+β3) -(β2+β4)，得到：

式(3)中dδasum,az為δasum,az的全微分函數，是垂向運動補償的殘余誤差，δaz和δβsum分別表示載體垂向加速度az和綜合安裝誤差角βsum的估計誤差。由于重力梯度信號通過機械旋轉的方式調制在敏感器測量信號的2ω頻率分量上，加速度計在2ω頻率分量上的測量精度約為水平，通過對加速度計敏感軸精密裝調，可將加速度計敏感軸垂向角度綜合安裝誤差控制在0.5′以內，即由此可計算載體垂向加速度測量誤差產生的重力梯度垂向運動補償殘余誤差在0.5E 水平，可忽略不計。

某次重力梯度船載測量實驗中，載體垂向加速度的時域圖和頻域圖如圖4 所示，其在2ω頻率分量上的幅值約為量級，以敏感器旋轉半徑100 mm 和濾波時間100 s 計算，為使垂向運動補償的殘余誤差控制在10E 水平，需要加速度計敏感軸垂向角度綜合安裝誤差的估計精度δβsum控制在0.5″以內，這對于高精度重力梯度動態測量造成了嚴峻挑戰。

圖4 船載實驗中載體垂向加速度的時域圖和頻域圖Fig.4 Time-domain and frequency-domain plots of carrier vertical acceleration in shipboard experiments

2 綜合安裝誤差角實時估計方法

綜合安裝誤差的估計精度δβsum要求控制在0.5″以內，只能利用載體垂向運動與敏感器輸出之間的相關性，實時估計該綜合安裝誤差角。

2.1 回歸方程構建

在式(1)的測量方程中考慮載體垂向加速度的影響，得：

通過遞推最小二乘的方法求取加速度計敏感軸垂向角度綜合安裝誤差，其遞推公式為：

式中，K是修正系數矩陣，P是狀態估計的均方誤差陣，是狀態的最優估計，I 是單位矩陣，R是正定的量測均方誤差陣，即遞推中的初值為：

利用載體垂向運動產生敏感器測量誤差的解析關系建立回歸方程，并通過遞推最小二乘的方法實時估計綜合安裝誤差角，從而結合載體垂向加速度完成重力梯度垂向運動補償，提高重力梯度實時測量精度。但遞推最小二乘方法是綜合利用了歷史所有量測值Zk，在經過較長時間遞推后，Pk陣逐漸收斂，修正系數矩陣Kk逐漸減小，使得新的量測值對狀態估計的修正作用逐漸減弱，致使狀態的最優估計不能很好地跟蹤狀態實際變化。此外，量測均方誤差陣R則與載體動態密切相關，同一數值無法滿足全測量周期內不同的載體運動擾動下濾波要求。而R陣的設置誤差會引起濾波精度降低，嚴重時還可導致濾波器發散，在實際工程應用中必須予以考慮。

2.2 含有遺忘因子與漸消因子的估計方法

針對上述問題，首先在遞推公式中引入遺忘因子，在狀態估計中刻意改變了遞推誤差與量測誤差之間的權重，從而提高了狀態估計對實際狀態的跟蹤能力?？蓪⑹?6)改寫為：

式(8)中s是略小于1 的實數比例因子，稱為遺忘因子。工程實現中可通過加速度計敏感軸垂向角度綜合安裝誤差的短期變化規律確定遺忘因子s的數值。

此外，為使得狀態估計遞推方程具備對不同工況的適應能力，在濾波中增加Sage-Husa 自適應環節[10]。量測矩陣的預測均方誤差矩陣可定義為：

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

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

式中初值λ0為1，c（0＜c＜ 1）是漸消因子，當k→∞時，有λk→ 1-c，此時可使濾波器始終擁有自適應能力。此外，通過序貫濾波的方法限制R陣的對角線元素數值，以保障R陣的正定性，滿足濾波器的實現要求。記則：

圖5 含有遺忘因子的Sage-Husa 自適應濾波算法流程圖Fig.5 Flowchart of Sage-Husa adaptive filtering algorithm with forgetting factor

3 遺忘因子和漸消因子的數值確定

3.1 遺忘因子的數值確定

在遞推最小二乘算法中引入遺忘因子s，能夠改善狀態估計的跟蹤性能，但也以s的比例提高了新量測信息的精度，等效降低了狀態估計算法的穩定性。因此，在遺忘因子s的參數選擇中，既要考慮狀態估計的跟蹤性能，又要考慮算法的魯棒性。本文通過半物理仿真的方法確定該參數的數值，其中綜合安裝誤差角βsum在動態下的測量規律由靜態下獲取，并結合某次船載重力梯度測量中的實測運動數據，構建仿真數據集。仿真中選取船載實驗實測數據海況為3 級，艦船排水量為3000 t，重力梯度測量過程中航速為10 knot。無需考慮綜合安裝誤差角實時估計中R陣的自適應問題，因此在本次仿真中將R陣設置為固定值。對該敏感器仿真信號使用遞推最小二乘法（s=1）以及含有遺忘因子典型值的遞推最小二乘方法對綜合安裝誤差角βsum進行狀態估計，多次仿真后某次典型結果如圖6 所示，待狀態估計收斂后，通過計算均方根（RMS）的方法評價狀態估計精度，相關統計結果如圖7 所示。結合圖6 和圖7 可知，當遺忘因子s取0.995 時，對綜合安裝誤差角βsum的估計精度最高，優于0.07″，因此在本研究中遺忘因子s取值為0.995。

圖6 不同遺忘因子情況下綜合安裝誤差角估計結果Fig.6 Estimation results of integrated installation error angle for different forgetting factor cases

3.2 漸消因子的數值確定

與遺忘因子的取值原則類似，漸消因子c的取值越小，對最新量測噪聲變化的適應能力越強，但若漸消因子c的取值過小，噪聲均方誤差陣的估計結果會劇烈變化，從而導致狀態估計精度降低，嚴重時會引起系統發散。一般來講，漸消因子c的取值范圍在0.9和0.999 之間。本文同樣采用半物理仿真的方法確定漸消因子的數值，仿真數據的構建方法與遺忘因子仿真數據的構建方法一致，但為有針對性提高量測噪聲的自適應能力，仿真中載體運動選取船載實驗實測數據海況為3-5 級，其中兩端數據部分海況為3 級，中間數據部分海況為5 級，遺忘因子s設置為0.995。

對該敏感器仿真信號使用不同漸消因子的最小二乘方法對綜合安裝誤差角βsum進行狀態估計，多次仿真后某次典型結果如圖8 所示，待狀態估計收斂后，通過計算均方根（RMS）的方法評價狀態估計精度，相關統計結果如圖9 所示。結合圖8 和圖9 可知，當漸消因子c取0.9385 時，對綜合安裝誤差角βsum的估計精度最高，優于0.23″，因此在本研究中漸消因子c的數值取0.9385。

圖8 不同漸消因子情況下綜合安裝誤差角估計結果Fig.8 Estimation results of the integrated installation error angle for different asymptotic factors

圖9 不同漸消因子情況下綜合安裝誤差角估計誤差Fig.9 Combined installation error angle estimation error with different asymptotic factors

4 實驗驗證

4.1 綜合安裝誤差估計結果

項目組在2022 年在南海某海域開展船載水平張量重力梯度測量實驗，實驗現場如圖10 所示。實驗中有效重力梯度測量時間為8 天，其中第3-4 天海況為5 級，其余時間海況為3-4 級。分別使用遞推最小二乘方法和本文提出的含有遺忘因子與自適應濾波的估計方法處理儀器實測數據，兩種方法得到的綜合安裝誤差角βsum估計值如圖11 所示。

圖10 船載水平張量重力梯度測量實驗現場圖Fig.10 Site map of shipboard horizontal tensor gravity gradient measurement experiment

圖11 兩種估計方法得到的 βsum估計結果Fig.11 βsumestimates obtained by both estimation methods

由圖11 可知，采用遞推最小二乘方法得到的βsum估計值無法反映βsum的時變特性，而采用遺忘因子與自適應濾波的估計方法得到的βsum估計值存在24 h的周期特性，能夠反映βsum隨溫度的變化特性，表明該方法能夠提高重力梯度儀對βsum的實時估計精度。

4.2 重力梯度測量結果及精度統計

分別使用遞推最小二乘方法和本文提出的含有遺忘因子與自適應濾波的估計方法處理儀器實測數據，得到兩路重力梯度測量信號重復線測量結果如圖12和圖13 所示。

圖12 遞推最小二乘方法下水平張量重力梯度測量結果Fig.12 Measurement results of horizontal tensor gravity gradient using recursive least squares method

圖13 遺忘因子與自適應濾波的估計方法下兩路重力梯度信號測量結果Fig.13 Measurements of two gravity gradient signals under the estimation method of forgetting factor and adaptive filtering

重力梯度動態測量通過內符合中誤差的數值來評價測量精度，其計算公式為[12]：

式(16)中，εj是重力梯度j分量的內符合中誤差，Γj1(i)和Γj2是兩次重力梯度動態測量j分量在第i點的測量值，n是測量總數。不同處理方法重力梯度測量重復線內符合精度統計結果如表1 所示，結果表明本文提出的估計方法能夠將內符合精度由30E@1km提高至15E@1km。

表1 不同補償方法重力梯度測量內符合精度統計Tab.1 Statistics on the internal compliance accuracy of gravity gradient measurements with different compensation methods

5 結論

本文針對基于旋轉加速度計原理的水平張量重力梯度儀在動態測量中加速度計敏感軸垂向綜合安裝誤差角的時變特性，以及狀態估計中易受到載體運動干擾的問題，提出了一種含有遺忘因子與自適應濾波的綜合安裝誤差實時估計方法。通過在狀態估計中增加遺忘因子的方式抑制綜合安裝誤差角的時變特性，并在濾波中引入漸消因子提高該方法對載體動態的適應性，實現了對綜合安裝誤差角進行次優估計，以此完成重力梯度信號的實時垂向運動補償。船載重力梯度測量實驗結果表明，與遞推最小二乘方法相比，本文所提方法實時得到的綜合安裝誤差角能夠反映其隨溫度的變化特性，可將重力梯度船載測量內符合精度由30E@1km 提高至15E@1km，具有工程指導意義。