一種用于重力梯度動態測量的載體環境引力梯度補償方法

嚴 飛，李 達，李 中,3，孟兆海,3,4,5，王懷君

（1.海裝北京局駐天津地區某軍事代表室，天津 300131;2.天津航海儀器研究所，天津 300131;3.中船集團航海保障技術實驗室，天津 300131;4 中國科學院地質與地球物理研究所，北京 100029;5 中國科學院油氣資源研究重點實驗室，北京 100029）

重力梯度是重力矢量的空間變化率，重力梯度儀是用于測量重力梯度的精密設備，基于Bell Aerospace公司提出的旋轉加速度計測量原理的重力梯度儀是迄今唯一實用的近地表動態重力梯度儀[1-3]。重力梯度儀的核心部件是重力梯度敏感器，動態測量時，重力梯度敏感器安裝在慣性穩定平臺上，穩定在地理坐標系下，測量地理坐標系下的重力梯度。在測量過程中，飛機、艦船等載體的姿態變化，將引起載體質量分布與重力梯度敏感器的相對位置發生變化。載體質量對重力梯度敏感器產生的引力梯度將隨著載體姿態變化而產生變化，形成載體環境引力梯度變化[4,5]。隨著重力梯度儀精度水平不斷提高，在動態測量中載體環境引力梯度變化已成為高精度重力梯度儀的一項重要誤差來源[6]。

本文基于重力梯度張量在不同坐標系下的轉換關系建立了載體環境引力梯度變化的解析模型，以此為基礎建立了回歸方程。通過Tikhonov正則化方法抑制了該方程解的不適定問題，得到載體環境引力梯度在載體坐標系下的數值解，以此進行載體環境引力梯度變化改正，通過實測海試數據驗證該方法的有效性，為動態重力梯度儀研制提供參考。

1 載體環境引力梯度變化模型

1.1 重力梯度儀測量原理

重力梯度是重力位的空間二階導數，表征重力矢量的空間變化率。在地理坐標系中，重力矢量可以分解為x，y，z三個方向上的三個分量中，每一分量沿平行于坐標軸方向均有一個梯度。因此，重力梯度張量共有9個分量,可將其表示為：

式中，Γij（i,j=x,y,z）為梯度張量的分量，表示重力分量gi在j方向上的變化率。

旋轉加速度計式重力梯度儀主要由重力梯度敏感器和慣性穩定平臺兩個關鍵組件構成。重力梯度敏感器主要用于完成重力梯度張量水平分量的測量；慣性穩定平臺則用于承載重力梯度敏感器，為重力梯度動態測量提供穩定的動力學環境。如圖1所示，重力梯度敏感器基于加速度計位置差分的測量原理，通過機械旋轉的方式將旋轉中心處的重力梯度張量水平分量調制到系統旋轉頻率的二倍頻處，加速度計四路和信號與重力梯度張量水平分量之間的關系可表示為[7]：

式中，a1、a2、a3、a4是四只加速度計敏感軸方向的加速度，l是加速度計檢測質心到旋轉中心的距離，Γxx、Γyy、Γxy是測量平面內對應方向的重力梯度張量分量，單位是 E（1 E=10-9s-2），ω是旋轉圓盤的旋轉角速率。進行動態測量時，慣性穩定平臺采用三環固定指北半解析式控制方案，在隔離載體角運動的同時，將重力梯度敏感器穩定在地理坐標系，為測量提供坐標基準。

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

1.2 重力梯度坐標變換

依據重力勢函數梯度和方向導數關系，得到重力梯度張量在不同坐標系下的變換公式為[8]：

式中：Γn是重力梯度張量在地理坐標系下的投影；Γb是重力梯度張量在載體坐標系下的投影；Cbn是從載體坐標系到地理坐標系的坐標變換矩陣。定義地理坐標系為“東-北-天”地理坐標系時，依據歐拉轉動原理，Cbn可表示為：

式中，ψ、θ和γ分別為載體的航向角、俯仰角和橫滾角。

重力場是一個保守場，表明重力場是一個無源場，且場線在空間內不閉合，有[9]：

式（5）表示重力場的旋度為零，說明重力梯度張量矩陣具有對稱性，即：

式（6）表示重力場的散度為零，意味著重力梯度張量矩陣的跡為零，即：

因此，得到重力梯度的坐標變換公式為：

1.3 載體環境引力梯度變化模型

式（10）可簡寫為：

式中，Rbn為從載體坐標系到地理坐標系的重力梯度張量列向量的坐標變換矩陣。

在動態測量過程中重力梯度敏感器始終穩定在地理坐標系，載體質量的分布與敏感器的相對位置發生變化，形成載體環境引力梯度變化。結合式（11），得到重力梯度儀載體環境引力梯度變化的解析方程為：

式中，δTen是由載體環境引力梯度變化引起的水平分量重力梯度測量誤差；Teb是載體質量在載體坐標系下對重力梯度敏感器旋轉中心位置產生的引力梯度。

利用實測艦船姿態信息結合艦船結構引力梯度正演計算，滿載排水量為5000 t的艦船在測線上載體環境引力梯度變化量約為30 E量級。

2 基于Tikhonov正則化的載體環境梯度誤差補償方法

2.1 回歸方程

依據式（12），可以寫出考慮載體環境引力梯度變化的重力梯度儀測量方程為：

利用載體環境引力梯度變化和載體姿態的解析關系建立回歸方程，可用多元回歸的方式求取載體質量在載體坐標系下對重力梯度敏感器旋轉中心位置引起的引力梯度變化，從而結合載體姿態信息對測量數據進行載體環境引力梯度變化改正，提高測量精度。

2.2 解的不適定分析

利用最小二乘法求解式（14），可得[11]：

式中是x的估計解，即載體質量在載體坐標系下對重力梯度敏感器旋轉中心位置引起引力梯度的估計值。

在不考慮回歸算子K的擾動時，有估計解的擾動式為[12]：

式中cond（K）是算子K在2范數意義下的條件數，其計算式為：

式（16）表明了回歸方程式（14）中的非線性量δ以cond（K）的倍數進行放大，因此為抑制擾動量δ對估計解的擾動，要求算子K的條件數越小越好。

在進行重力梯度測量作業時，通常由飛機或艦船搭載重力梯度儀，在固定區域內沿網格化測線測量。在測線中保持勻速直線運動，最大程度降低載體加速度對于重力梯度測量的干擾。因此在測量過程中載體俯仰角和橫滾角的波動較小，只有航向角能產生較大變化。而在飛機爬升和艦船轉彎等俯仰角、橫滾角波動較大的航行狀態下，載體運動加速度干擾同樣很大，此時重力梯度儀無法實現高精度測量。

某次重力梯度儀船載試驗中兩條反向的測線船體俯仰角、橫滾角和航向角如圖2，3所示（定義為測線1#和測線2#），兩條測線的航跡線如圖4所示。從中可以看出測量過程中船體俯仰角波動峰峰值為2°，橫滾角波動峰峰值為2°，航向角波動峰峰值為5°。計算cond（K）的結果如表1所示。

圖2 重力梯度儀船載試驗中測線1#船體航姿曲線Fig.2 Line 1# ship attitude curve in the shipboard test of gravity gradiometer

圖3 重力梯度儀船載試驗中測線2#船體航姿曲線Fig.3 Line 2# ship attitude curve in the shipboard test of gravity gradiometer

圖4 重力梯度儀船載試驗中測線1#和測線2#航跡線Fig.4 The track lines of line 1# and line 2# in the shipborne test of gravity gradiometer

表1 測線1#和測線2#的cond(K)計算結果Tab.1 Calculation results of cond(K)of line 1# and line 2# respectively

從表1中可以看出，利用兩條反向的測線1#和2#上的船體姿態數據，cond（K）為3×103～5×103量級，仍是病態的，即具有嚴重的復共線性特性。因此考慮利用兩條正交的測線進行回歸分析，降低cond（K）的數值。

引入一條與測線1#正交的測線，定義為測線3#，其船體俯仰角、橫滾角和航向角如圖5所示，測線1#和測線3#的航跡線如圖6所示，利用1#、2#和3#測線的數據計算cond(K)的結果如表2所示。

圖5 重力梯度儀船載試驗中測線3#船體航姿曲線Fig.5 Line 3# ship attitude curve in the shipboard test of gravity gradiometer

圖6 重力梯度儀船載試驗中測線1#和測線3#航跡線Fig.6 The track lines of line 1# and line 3# in the shipborne test of gravity gradiometer

表2 測線1#和測線3#的cond(K)計算結果Tab.2 Calculation results of cond(K)of line 1# and line 3# respectively

從表2中可以看出，即使利用3條正交的測線（1#、2#和3#）上的船體姿態數據計算cond（K），數值仍在1×103量級，仍是病態的，在式（14）中會放大測量信號yδ中的擾動量δ。同時，回歸方程的擾動量δ不僅包含測量噪聲，還包含真實的地表重力梯度信號，回歸方程的擾動量是擬合量的5倍以上，為防止“過擬合”的情況，必須選擇合理的正則化方法，以求取回歸方程的精確估計解。

2.3 基于Tikhonov正則化的回歸方法

目前解決最小二乘回歸病態問題應用較廣的方法是Tikhonov正則化方法[13]，它通過求解如式（18）所示的最小化問題獲得正則數值解。

可通過式（19）得到。

式中，I是單位矩陣，α是控制參數（α＞ 0）。其中（KTK+αI ）-1KT是正則化算子。Tikhonov正則化方法的基本思想是犧牲估計解的無偏性將系數矩陣K的病態問題轉化為良性問題，使得正則數值解更加逼近于精確解x。

2.4 正則化參數確定

目前已有大量文獻研究了Tikhonov正則化方法中正則化算子的構造和相關正則化參數的選取[14,15]，本文采用半物理仿真的方式選擇控制參數α的最佳數值。通過引力梯度正演的方法仿真得到該方程的理論解x，利用實測重力梯度數據和實測姿態信息結合式（14）分別計算δ、K和yδ。定義誤差比為?？刂茀郸猎冢?.0001,100］范圍內分別求取回歸方程的估計解，計算步長為0.0001，并計算誤差比，結果如圖7所示。作為對比，再使用常規最小二乘的方法求取估計解，并計算誤差比。典型計算結果見表3。

圖7 不同控制參數α時的誤差比Fig.7 Error ratio of different control parameter α

表3 不同條件下回歸方程估計解的計算結果Tab.3 Calculation results of estimated solution of regression equation under different conditions

表3 不同條件下回歸方程估計解的計算結果Tab.3 Calculation results of estimated solution of regression equation under different conditions

估計值/E b Γ b xx Γ b xy Γ b xz Γ b yy Γ 誤差比/%yz理論值 26.2225.874.01 -5.332.88 / α估=計解 860.7225.563.04829.86 -50.7823.220.0001估α計=解 15.9225.443.50 -13.909.681.989 α估計=解 15.3725.280.94 -14.594.331.9926 α估=計解 11.9422.27-0.71 -11.900.6515.141000最估小計二解乘 878.1725.552.71847.26 -52.6924.69

結合圖7和表3中可知，與較常規最小二乘方法比，Tikhonov正則化的方法只要將控制參數α選取在9～26之間，能將補償量的誤差比從24.69%降低至2%以內，且相差不大。因此在本研究中選擇α= 15。

3 試驗驗證

3.1 載體環境引力梯度估計

重力梯度儀船載試驗中測線1# 和2# 兩條重復線的測量數據未進行載體環境引力梯度補償時，重力梯度測量信號結果如圖8和9所示。結合載體姿態數據使用式（19）進行載體環境梯度估計，分別使用常規最小二乘法和Tikhonov正則化方法進行回歸估計，計算結果如表4所示。

圖8 重力梯度 Γuv分量原始輸出信號Fig.8 Original output signal of gravity gradient component Γuv

圖9 重力梯度 Γxy分量原始輸出信號Fig.9 Original output signal of gravity gradient component Γxy

表4 載體環境引力梯度估計結果Tab.4 Calculation results of carrier environmental gradient

3.1 重力梯度測量結果及精度評價

將Tikhonov正則化法和常規最小二乘法得到的載體環境梯度估計值代入式（14），進行重力梯度信號載體環境引力梯度變化補償，使用Tikhonov正則化法補償后重力梯度測量信號結果如圖10和11所示，使用最小二乘法補償后重力梯度測量信號結果如圖12和13所示，分別計算補償前后重力梯度信號內符合精度，如表5所示。

圖10 Tikhonov正則化法載體環境引力梯度變化補償后重力梯度分量信號Fig.10. component signal of gravity gradient after compensation of environmental gravity gradient change by Tikhonov regularization method

圖11 Tikhonov正則化法載體環境引力梯度變化補償后重力梯度 Γxy分量信號Fig.11 Γxy component signal of gravity gradient after compensation of environmental gravity gradient change by Tikhonov regularization method

圖12 最小二乘法載體環境引力梯度變化補償后 重力梯度 Γxy分量信號Fig.12 Γuv component signal of gravity gradient after compensation of gravity gradient change of carrier environment by least square method

圖13 最小二乘法法載體環境引力梯度變化補償后 重力梯度 Γxy分量信號Fig.13.Γxy component signal of gravity gradient after compensation of gravity gradient change of carrier environment by least square method

表5 載體環境引力梯度變化補償前后內符合中誤差對比Tab.5 Comparison of the mean square error of the internal coincidence before and after compensation for the change of the gravity gradient of the carrier environment

4 結 論

本文針對重力梯度儀動態測量過程中載體環境引力梯度干擾的問題，建立了動態測量過程中載體環境引力梯度變化的傳遞模型，以此為基礎建立了回歸方程。通過Tikhonov正則化方法抑制了該方程解的不適定問題，利用半實物仿真的方法得到最優控制參數α的取值范圍，使補償量的誤差控制在2%。經實測船載數據驗證，該方法切實可行，補償后內符合中誤差較原始重力梯度測量信號中誤差分別降低19 E和 21 E，較常規最小二乘法進一步降低7 E和3 E，補償后重力梯度儀兩路測量信號內符合中誤差均達到優于10 E的精度水平。