一種用于重力梯度動(dòng)態(tài)測(cè)量的載體環(huán)境引力梯度補(bǔ)償方法

嚴(yán) 飛，李 達(dá)，李 中,3，孟兆海,3,4,5，王懷君

（1.海裝北京局駐天津地區(qū)某軍事代表室，天津 300131;2.天津航海儀器研究所，天津 300131;3.中船集團(tuán)航海保障技術(shù)實(shí)驗(yàn)室，天津 300131;4 中國(guó)科學(xué)院地質(zhì)與地球物理研究所，北京 100029;5 中國(guó)科學(xué)院油氣資源研究重點(diǎn)實(shí)驗(yàn)室，北京 100029）

重力梯度是重力矢量的空間變化率，重力梯度儀是用于測(cè)量重力梯度的精密設(shè)備，基于Bell Aerospace公司提出的旋轉(zhuǎn)加速度計(jì)測(cè)量原理的重力梯度儀是迄今唯一實(shí)用的近地表動(dòng)態(tài)重力梯度儀[1-3]。重力梯度儀的核心部件是重力梯度敏感器，動(dòng)態(tài)測(cè)量時(shí)，重力梯度敏感器安裝在慣性穩(wěn)定平臺(tái)上，穩(wěn)定在地理坐標(biāo)系下，測(cè)量地理坐標(biāo)系下的重力梯度。在測(cè)量過(guò)程中，飛機(jī)、艦船等載體的姿態(tài)變化，將引起載體質(zhì)量分布與重力梯度敏感器的相對(duì)位置發(fā)生變化。載體質(zhì)量對(duì)重力梯度敏感器產(chǎn)生的引力梯度將隨著載體姿態(tài)變化而產(chǎn)生變化，形成載體環(huán)境引力梯度變化[4,5]。隨著重力梯度儀精度水平不斷提高，在動(dòng)態(tài)測(cè)量中載體環(huán)境引力梯度變化已成為高精度重力梯度儀的一項(xiàng)重要誤差來(lái)源[6]。

本文基于重力梯度張量在不同坐標(biāo)系下的轉(zhuǎn)換關(guān)系建立了載體環(huán)境引力梯度變化的解析模型，以此為基礎(chǔ)建立了回歸方程。通過(guò)Tikhonov正則化方法抑制了該方程解的不適定問(wèn)題，得到載體環(huán)境引力梯度在載體坐標(biāo)系下的數(shù)值解，以此進(jìn)行載體環(huán)境引力梯度變化改正，通過(guò)實(shí)測(cè)海試數(shù)據(jù)驗(yàn)證該方法的有效性，為動(dòng)態(tài)重力梯度儀研制提供參考。

1 載體環(huán)境引力梯度變化模型

1.1 重力梯度儀測(cè)量原理

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

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

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

式中，a1、a2、a3、a4是四只加速度計(jì)敏感軸方向的加速度，l是加速度計(jì)檢測(cè)質(zhì)心到旋轉(zhuǎn)中心的距離，Γxx、Γyy、Γxy是測(cè)量平面內(nèi)對(duì)應(yīng)方向的重力梯度張量分量，單位是 E（1 E=10-9s-2），ω是旋轉(zhuǎn)圓盤的旋轉(zhuǎn)角速率。進(jìn)行動(dòng)態(tài)測(cè)量時(shí)，慣性穩(wěn)定平臺(tái)采用三環(huán)固定指北半解析式控制方案，在隔離載體角運(yùn)動(dòng)的同時(shí)，將重力梯度敏感器穩(wěn)定在地理坐標(biāo)系，為測(cè)量提供坐標(biāo)基準(zhǔn)。

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

1.2 重力梯度坐標(biāo)變換

依據(jù)重力勢(shì)函數(shù)梯度和方向?qū)?shù)關(guān)系，得到重力梯度張量在不同坐標(biāo)系下的變換公式為[8]：

式中：Γn是重力梯度張量在地理坐標(biāo)系下的投影；Γb是重力梯度張量在載體坐標(biāo)系下的投影；Cbn是從載體坐標(biāo)系到地理坐標(biāo)系的坐標(biāo)變換矩陣。定義地理坐標(biāo)系為“東-北-天”地理坐標(biāo)系時(shí)，依據(jù)歐拉轉(zhuǎn)動(dòng)原理，Cbn可表示為：

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

重力場(chǎng)是一個(gè)保守場(chǎng)，表明重力場(chǎng)是一個(gè)無(wú)源場(chǎng)，且場(chǎng)線在空間內(nèi)不閉合，有[9]：

式（5）表示重力場(chǎng)的旋度為零，說(shuō)明重力梯度張量矩陣具有對(duì)稱性，即：

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

因此，得到重力梯度的坐標(biāo)變換公式為：

1.3 載體環(huán)境引力梯度變化模型

式（10）可簡(jiǎn)寫為：

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

在動(dòng)態(tài)測(cè)量過(guò)程中重力梯度敏感器始終穩(wěn)定在地理坐標(biāo)系，載體質(zhì)量的分布與敏感器的相對(duì)位置發(fā)生變化，形成載體環(huán)境引力梯度變化。結(jié)合式（11），得到重力梯度儀載體環(huán)境引力梯度變化的解析方程為：

式中，δTen是由載體環(huán)境引力梯度變化引起的水平分量重力梯度測(cè)量誤差；Teb是載體質(zhì)量在載體坐標(biāo)系下對(duì)重力梯度敏感器旋轉(zhuǎn)中心位置產(chǎn)生的引力梯度。

利用實(shí)測(cè)艦船姿態(tài)信息結(jié)合艦船結(jié)構(gòu)引力梯度正演計(jì)算，滿載排水量為5000 t的艦船在測(cè)線上載體環(huán)境引力梯度變化量約為30 E量級(jí)。

2 基于Tikhonov正則化的載體環(huán)境梯度誤差補(bǔ)償方法

2.1 回歸方程

依據(jù)式（12），可以寫出考慮載體環(huán)境引力梯度變化的重力梯度儀測(cè)量方程為：

利用載體環(huán)境引力梯度變化和載體姿態(tài)的解析關(guān)系建立回歸方程，可用多元回歸的方式求取載體質(zhì)量在載體坐標(biāo)系下對(duì)重力梯度敏感器旋轉(zhuǎn)中心位置引起的引力梯度變化，從而結(jié)合載體姿態(tài)信息對(duì)測(cè)量數(shù)據(jù)進(jìn)行載體環(huán)境引力梯度變化改正，提高測(cè)量精度。

2.2 解的不適定分析

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

式中是x的估計(jì)解，即載體質(zhì)量在載體坐標(biāo)系下對(duì)重力梯度敏感器旋轉(zhuǎn)中心位置引起引力梯度的估計(jì)值。

在不考慮回歸算子K的擾動(dòng)時(shí)，有估計(jì)解的擾動(dòng)式為[12]：

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

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

在進(jìn)行重力梯度測(cè)量作業(yè)時(shí)，通常由飛機(jī)或艦船搭載重力梯度儀，在固定區(qū)域內(nèi)沿網(wǎng)格化測(cè)線測(cè)量。在測(cè)線中保持勻速直線運(yùn)動(dòng)，最大程度降低載體加速度對(duì)于重力梯度測(cè)量的干擾。因此在測(cè)量過(guò)程中載體俯仰角和橫滾角的波動(dòng)較小，只有航向角能產(chǎn)生較大變化。而在飛機(jī)爬升和艦船轉(zhuǎn)彎等俯仰角、橫滾角波動(dòng)較大的航行狀態(tài)下，載體運(yùn)動(dòng)加速度干擾同樣很大，此時(shí)重力梯度儀無(wú)法實(shí)現(xiàn)高精度測(cè)量。

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

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

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

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

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

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

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

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

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

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

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

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

目前解決最小二乘回歸病態(tài)問(wèn)題應(yīng)用較廣的方法是Tikhonov正則化方法[13]，它通過(guò)求解如式（18）所示的最小化問(wèn)題獲得正則數(shù)值解。

可通過(guò)式（19）得到。

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

2.4 正則化參數(shù)確定

目前已有大量文獻(xiàn)研究了Tikhonov正則化方法中正則化算子的構(gòu)造和相關(guān)正則化參數(shù)的選取[14,15]，本文采用半物理仿真的方式選擇控制參數(shù)α的最佳數(shù)值。通過(guò)引力梯度正演的方法仿真得到該方程的理論解x，利用實(shí)測(cè)重力梯度數(shù)據(jù)和實(shí)測(cè)姿態(tài)信息結(jié)合式（14）分別計(jì)算δ、K和yδ。定義誤差比為?？刂茀?shù)α在［0.0001,100］范圍內(nèi)分別求取回歸方程的估計(jì)解，計(jì)算步長(zhǎng)為0.0001，并計(jì)算誤差比，結(jié)果如圖7所示。作為對(duì)比，再使用常規(guī)最小二乘的方法求取估計(jì)解，并計(jì)算誤差比。典型計(jì)算結(jié)果見(jiàn)表3。

圖7 不同控制參數(shù)α?xí)r的誤差比Fig.7 Error ratio of different control parameter α

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

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

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

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

3 試驗(yàn)驗(yàn)證

3.1 載體環(huán)境引力梯度估計(jì)

重力梯度儀船載試驗(yàn)中測(cè)線1# 和2# 兩條重復(fù)線的測(cè)量數(shù)據(jù)未進(jìn)行載體環(huán)境引力梯度補(bǔ)償時(shí)，重力梯度測(cè)量信號(hào)結(jié)果如圖8和9所示。結(jié)合載體姿態(tài)數(shù)據(jù)使用式（19）進(jìn)行載體環(huán)境梯度估計(jì)，分別使用常規(guī)最小二乘法和Tikhonov正則化方法進(jìn)行回歸估計(jì)，計(jì)算結(jié)果如表4所示。

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

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

表4 載體環(huán)境引力梯度估計(jì)結(jié)果Tab.4 Calculation results of carrier environmental gradient

3.1 重力梯度測(cè)量結(jié)果及精度評(píng)價(jià)

將Tikhonov正則化法和常規(guī)最小二乘法得到的載體環(huán)境梯度估計(jì)值代入式（14），進(jìn)行重力梯度信號(hào)載體環(huán)境引力梯度變化補(bǔ)償，使用Tikhonov正則化法補(bǔ)償后重力梯度測(cè)量信號(hào)結(jié)果如圖10和11所示，使用最小二乘法補(bǔ)償后重力梯度測(cè)量信號(hào)結(jié)果如圖12和13所示，分別計(jì)算補(bǔ)償前后重力梯度信號(hào)內(nèi)符合精度，如表5所示。

圖10 Tikhonov正則化法載體環(huán)境引力梯度變化補(bǔ)償后重力梯度分量信號(hào)Fig.10. component signal of gravity gradient after compensation of environmental gravity gradient change by Tikhonov regularization method

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

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

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

表5 載體環(huán)境引力梯度變化補(bǔ)償前后內(nèi)符合中誤差對(duì)比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 結(jié) 論

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