基于球諧模型與多傳感器融合的高精度重力擾動補償方法

劉宇鑫 王新龍 王勛 高文寧 胡曉東

引用格式：劉宇鑫，王新龍，王勛，等.基于球諧模型與多傳感器融合的高精度重力擾動補償方法［J］.航空兵器，2023，30（1）：104-113.

LiuYuxin，WangXinlong，WangXun，etal.AnAccurateGravityDisturbanceCompensationMethodBasedonSphericalHarmonicModelandMultiSensorFusion［J］.AeroWeaponry，2023，30（1）：104-113.（inChinese）

摘要：隨著對高精度長航時慣導系統(tǒng)性能要求的提升，重力擾動已成為慣導系統(tǒng)的主要誤差源，能否對其進行有效補償是提升導航精度的關(guān)鍵因素。傳統(tǒng)基于球諧模型的補償方法無法反映地球重力場的細節(jié)信息，對中短波重力擾動分量補償效果不佳；傳統(tǒng)狀態(tài)估計法中馬爾科夫狀態(tài)模型對長波重力擾動分量描述精度較差，因此對長波分量的補償效果不佳。以上方法均無法對波段較寬的實際重力擾動進行精確補償，針對這一問題，本文設(shè)計了一種綜合利用球諧模型與多傳感器信息融合的高精度重力擾動補償方法。該方法一方面利用低階球諧模型對長波重力擾動分量計算精度較高的特點，對長波分量進行補償；另一方面利用捷聯(lián)慣導/激光多普勒測速/氣壓高度計構(gòu)成完全自主的組合導航系統(tǒng)，將殘余的中短波重力擾動分量建立為高精度的馬爾科夫模型，從而實現(xiàn)對中短波分量的狀態(tài)估計與補償。仿真結(jié)果表明，所提出的重力擾動補償方法可有效改善重力擾動估計效果，進而實現(xiàn)高精度的重力擾動補償，具有較高的導航精度。

關(guān)鍵詞：捷聯(lián)慣導系統(tǒng)；重力擾動補償；重力場球諧模型；狀態(tài)估計；組合導航

中圖分類號：TJ765;V249.32+2

文獻標識碼：A

文章編號：1673-5048（2023）01-0104-10

DOI：10.12132/ISSN.1673-5048.2022.0038

0引言

慣性導航系統(tǒng)具有高度的自主性、隱蔽性以及信息完備等特點，目前廣泛應(yīng)用與軍事、工程、科學研究以及民用領(lǐng)域中［1］。在慣導解算過程中，通常采用理論重力模型計算理論重力并從比力測量值中對其進行補償［2］。但理論重力與實際重力并不相等，該差異即為重力擾動，其導致的補償殘差會造成慣導解算誤差。多數(shù)地區(qū)重力擾動的大小在15～80mGal之間［3］，是高精度慣導系統(tǒng)的主要誤差源之一［4-5］，對導航精度有較大的影響。

目前重力補償方法可根據(jù)是否已知載體運動區(qū)域的重力網(wǎng)格數(shù)據(jù)而分為兩大類。在已知區(qū)域重力網(wǎng)格數(shù)據(jù)庫時，通過對該數(shù)據(jù)庫進行插值獲取當前位置的重力擾動矢量并對其進行補償［6-8］。這種方法在局部區(qū)域內(nèi)的補償精度較高，且具有較好的實時性，但無法在該數(shù)據(jù)庫未覆蓋的區(qū)域適用。在工程中高精度的重力擾動測量數(shù)據(jù)往往較難獲取，因此這種補償方法的應(yīng)用區(qū)域受限。在缺乏區(qū)域重力網(wǎng)格數(shù)據(jù)庫時，重力擾動補償方法又可細分為兩種：基于重力場模型的補償方法以及基于狀態(tài)估計的補償方法。

在基于重力場模型的補償方法中，通常利用EGM2008球諧模型計算重力擾動并進行補償。這種方法適用于全球任意位置的重力擾動補償，但其計算結(jié)果無法反映地球重力場的細節(jié)信息［9-10］，因此導致該方法在山地等重力擾動變化較為顯著區(qū)域的補償精度較低，以中國南部高原區(qū)為例，其計算誤差的標準差達到20.74mGal［11］，難以滿足高精度的補償需求。此外，高階球諧模型的計算負擔較大，模型參數(shù)占用存儲空間大，計算耗時長［12］，不適用于傳統(tǒng)導航系統(tǒng)的硬件配置環(huán)境，并且難以滿足慣導系統(tǒng)解算實時性的要求。

狀態(tài)估計法將捷聯(lián)慣導系統(tǒng)與其他導航傳感器結(jié)合成組合導航系統(tǒng)，通過引入重力擾動狀態(tài)量，建立其狀態(tài)空間模型，利用最優(yōu)估計算法實現(xiàn)對重力擾動的估計與補償［13-14］。這類方法在補償重力擾動的同時，能有效抑制慣導系統(tǒng)誤差的發(fā)散，因此具有較好的補償效果。但在狀態(tài)估計過程中，重力擾動與慣導系統(tǒng)誤差耦合，需要借助高精度的重力擾動狀態(tài)模型才能對其進行有效估計，即重力擾動的建模精度是制約該方法重力補償精度的重要因素［15］。實際重力擾動可視為一系列不同幅值和波長的重力擾動分量信號的疊加，其所覆蓋的波段很廣，包含波長約1～10000km的分量信號［16］。然而常用的馬爾科夫重力擾動模型所能覆蓋的波段有限，當相關(guān)距離取幾十千米時，能對100km以下的重力擾動分量進行較為精確的描述［15，17］，但對長波分量的描述效果不佳，因此狀態(tài)估計法通常只能對中短波重力擾動分量進行較為有效的估計，對全波段的重力擾動估計精度有限。

針對以上問題，本文設(shè)計了一種球諧模型補償與狀態(tài)估計補償相互結(jié)合的重力擾動補償方法，在缺乏區(qū)域重力網(wǎng)格數(shù)據(jù)庫時能完全自主地實現(xiàn)重力擾動的估計與補償，并對該方法的性能進行了驗證。

1重力擾動對慣導解算影響分析

1.1重力擾動成因

在地球物理學中通常人為地選擇一個形狀規(guī)則、密度均勻的旋轉(zhuǎn)橢球體（即參考橢球）對地球做近似，由該參考橢球產(chǎn)生的重力為理論重力，記為γ0。以常用的1980年國際大地測量參考橢球為例，理論重力加速度計算公式為［18］

γ0=9.780325×［1+a1sin2L+a2sin2（2L）］－a3h（1）

式中：a1=0.00530240；a2=-0.00000582；a3=0.00003086；L為地理緯度；h為載體高度。

不過實際上地球并非是理想的旋轉(zhuǎn)橢球體，其形狀與參考橢球并不完全相同，而且地球內(nèi)部物質(zhì)分布結(jié)構(gòu)不規(guī)則，導致局部地區(qū)的密度不均，從而使得某一點處的實際重力加速度g與理論重力加速度γ0存在一定差異，由此可將這一差異定義為重力擾動δg，即

δg=g－γ0（2）

為了便于研究，通常將某點處實際重力方向與理論重力方向之間的夾角稱為垂線偏差，將其在卯酉圈上的分量記為η，在子午圈方向上的分量記為ξ；將重力擾動矢量在理論重力矢量方向上的分量稱為重力異常，記為δgU，如圖1所示。

由幾何關(guān)系可得重力擾動矢量在地理坐標系中的投影：

δg=［－γ0η－γ0ξ－δgU］T（3）

由此可見，重力擾動會引起實際重力與理論重力在當?shù)氐乩碜鴺讼迪氯齻€方向的誤差。

1.2重力擾動對慣導影響機理

捷聯(lián)慣導系統(tǒng)利用固聯(lián)于載體上的加速度計和陀螺儀分別測量載體相對慣性空間的比力f～b和角速度ω～bib，并以此為基礎(chǔ)進行導航解算，其原理如圖2所示。

在解算姿態(tài)時，利用陀螺測得角速度結(jié)合載體的位置和速度計算得到姿態(tài)角速度，進而通過姿態(tài)微分方程實時更新載體本體系b系與導航坐標系n系之間的坐標轉(zhuǎn)換矩陣C^nb，s（為與后續(xù)其他系統(tǒng)進行區(qū)分，利用下標s表示SINS相關(guān)計算參數(shù)），從而解算得到姿態(tài)角。在解算位置時，對加速度測量所得比力中的重力等有害加速度進行補償?shù)玫捷d體運動加速度，經(jīng)過兩次積分運算后可分別解算得到當前的速度和位置。

考慮慣性器件測量誤差以及初始誤差的影響，可得慣導誤差方程：

δv·ns=f^n×φns－（2ω^nie，s+ω^nen，s）×δvns－（2δωnie+δωnen）×v^ns－δgn0+C^nb，sδfb（4）

δL·s=δvN，sRM+h^s－v^N，s（RM+h^s）2δhs

δλ·s=secL^sRN+h^sδvE，s+v^E，ssecL^stanL^sRN+h^sδLs－v^E，ssecL^s（RN+h^s）2δhsδh·s=δvU，s（5）

φ·ns=δωnie，s－ω^nin，s×φns+δωnen，s－Cnb，sδωbib，s（6）

式中：RM為地球的子午圈半徑；RN為地球的卯酉圈半徑；ω^nie，s為地球自轉(zhuǎn)角速度；δωnie，s為地球自轉(zhuǎn)角速度的計算誤差；ω^nen，s為導航系相對地球系的轉(zhuǎn)動角速度；δωnie，s為導航系相對地球系的轉(zhuǎn)動角速度的計算誤差；δvns=［δvE，sδvN，sδvU，s］為速度誤差；δpns=［δLsδλsδhs］T為位置誤差；φs為姿態(tài)失準角；δg0為理論重力加速度的計算誤差；δfb為比力測量誤差；δωbib為陀螺儀的測量誤差。

然而由于重力擾動的影響，慣導工作時用于補償比力測量值中有害加速度的理論重力加速度與實際重力加速度并不相等，實際重力加速度為理論重力加速度與重力擾動之和，即

gn=gn0+δgn0（7）

僅采用理論重力加速度公式對比力測量值進行補償，會導致計算結(jié)果中存在一定的重力擾動殘差，進而引起額外的加速度解算誤差，此時加速度解算誤差δv～·ns為

δv～·ns=δv·ns－δgn（8）

將式（4）代入式（8），則可得考慮重力擾動時加速度解算誤差的具體形式：

δv～·ns=f^n×φns－（2ω^nie，s+ω^nen，s）×δvns－δgn－

C^nb，s

Δb－（2δωnie，s+δωnen，s）×vns－δgn0（9）

式中：Δb為加速度計測量誤差。

受重力擾動影響的加速度誤差經(jīng)過兩次積分后分別引起速度誤差δv^ns與位置誤差δp^ns。速度和位置誤差又進一步引起角速度計算誤差δω^nie和δω^nen：

δω^nie=0－ωiesinL^s·δL^s

ωiecosL^s·δL^s（10）

δω^nen=－δv^N，sRM+h^s+v^N，sδh^s（RM+h^s）2δv^E，sRN+h^s－v^E，sδh^s（RN+h^s）2

δv^E，stanL^sRN+h^s+v^E，sδL^ssec2L^sRN+h^s－v^E，sδh^stanL^s（RN+h^s）2（11）

由式（8）～（11）可知，重力擾動所引起額外的加速度解算誤差會導致額外的角速度計算誤差，進而通過式（6）引起額外的姿態(tài)解算誤差。姿態(tài)誤差又進一步通過式（4）～（5）造成位置誤差和速度誤差，形成閉環(huán)回路，持續(xù)對慣導系統(tǒng)的位置、速度和姿態(tài)解算產(chǎn)生影響。重力擾動對捷聯(lián)慣導系統(tǒng)的影響機理與加速度計零偏類似，其所造成的慣導解算誤差隨時間積累。由以上分析可知，重力擾動是影響高精度長航時慣導系統(tǒng)性能的關(guān)鍵因素之一。

2球諧模型與狀態(tài)估計相結(jié)合的補償方案設(shè)計

眾所周知，捷聯(lián)慣導系統(tǒng)通常采用理論重力加速度模型對比力測量值中的實際重力加速度進行粗略補償，因此其導航解算結(jié)果中包含由理論與實際重力差異（即重力擾動）所引起的導航誤差。

為了實現(xiàn)重力擾動的估計與補償，可借助其他導航傳感器提供不受重力擾動影響的導航參數(shù)，利用狀態(tài)估計算法從導航解算誤差中分離得到重力擾動估計值，并對其進行補償。與此同時，針對狀態(tài)估計法中馬爾科夫狀態(tài)模型對長波重力擾動描述精度有限的問題，考慮到導航系統(tǒng)的硬件配置限制，截取低階EGM2008球諧模型對長波的重力擾動進行補償，從而提升整體的補償效果。

2.1球諧模型補償

球諧模型補償部分根據(jù)EGM2008模型計算當前解算所得位置處的長波重力擾動分量，并在補償有害加速度時直接對加速度解算值進行校正，其原理如圖3所示。

通常用空間分辨率λres表示球諧模型對重力擾動描述的精細程度，其計算公式如下：

λres=πaN（12）

式中：a為地球的半長軸長度；N為球諧函數(shù)的最高階數(shù)。球諧模型的最高階數(shù)越大，相應(yīng)計算結(jié)果的分辨率越高，但受導航平臺硬件存儲空間以及計算實時性的限制，EGM2008模型階數(shù)小于12階時才能滿足導航系統(tǒng)計算資源要求［5］，因此截取前12階的參數(shù)計算重力擾動的長波分量，其具體計算公式可參見文獻［19］。

2.2狀態(tài)估計補償

通過將SINS與其他不受重力擾動影響的導航傳感器結(jié)合形成組合導航系統(tǒng)，對慣導系統(tǒng)誤差進行估計，同時通過最優(yōu)估計方法實現(xiàn)中短波重力擾動分量的估計與補償。

2.2.1陀螺/激光多普勒測速系統(tǒng)

在其他導航傳感器的選擇上，激光多普勒測速儀（LDV）具有測量精度高、動態(tài)響應(yīng)快、測量范圍大、測量線性度好等優(yōu)點，能提供高精度的速率測量信息［20］，通過與陀螺儀組合形成的導航系統(tǒng)可提供完備的導航信息，其原理如圖4所示。

后文中用下標D表示陀螺/多普勒測速儀系統(tǒng)的相關(guān)計算參數(shù)，以與SINS相關(guān)計算參數(shù)進行區(qū)分。利用激光多普勒測速儀測量載體相對地面的運動速率，利用陀螺儀測量所得角速度實時計算載體的姿態(tài)矩陣C^nb，D，再根據(jù)已知的安裝矩陣Cbm，將多普勒測速儀在其安裝坐標系下所測的速度v～mD投影到地理坐標系下得到

v^nD=v^E，Dv^N，Dv^U，DT=C^nb，DCbmv～mD（13）

由此可根據(jù)位置微分方程解算得到位置信息：

p^nD=L^Dλ^Dh^DT（14）

式中：L^·D=v^N，DRM+h^D；λ^·D=secL^DRN+h^Dv^E，D；h^·D=v^U，D。

與捷聯(lián)慣導系統(tǒng)類似，利用陀螺儀測量所得角速度ω～bib以及由位置和速度計算所得的ω^nie，D和ω^nen，D，可對姿態(tài)矩陣C^nb，D進行更新，即

C^·nb，D=C^nb，D（Ω～bib）－（Ω^nie，D+Ω^nen，D）C^nb，D（15）

式中：Ω～bib，Ω^nie，D和Ω^nen，D分別為ω～bib，ω^nie，D和ω^nen，D構(gòu)成的反對稱矩陣。

由式（13）～（15）可知，陀螺/激光多普勒測速導航系統(tǒng)直接利用激光多普勒測速儀測得的速度進行解算，不涉及重力的補償，因此其導航結(jié)果不受重力擾動的影響，可在狀態(tài)估計法中與慣導解算結(jié)果進行匹配，用于重力擾動的估計。

2.2.2狀態(tài)估計補償系統(tǒng)

利用SINS與陀螺/激光多普勒測速系統(tǒng)構(gòu)成組合導航系統(tǒng)，與此同時，由于兩導航系統(tǒng)均依賴陀螺的測量值，誤差會隨時間累積，因此用氣壓高度計對系統(tǒng)進行阻尼，其原理如圖5所示。

將SINS和陀螺/LDV系統(tǒng)解算所得導航參數(shù)之差以及氣壓高度計提供的高度信息作為量測值，同時建立中短波重力擾動分量的狀態(tài)模型，進而利用最優(yōu)估計算法對重力擾動以及導航誤差進行狀態(tài)估計與補償。

2.3重力擾動補償方案設(shè)計

綜合球諧模型以及狀態(tài)估計補償法，可得如圖6所示的重力擾動補償方案。

所設(shè)計方案同時采取球諧模型以及狀態(tài)估計兩種途徑對重力擾動進行補償。由解算所得位置，利用球諧模型計算長波重力擾動分量并進行補償；基于組合導航系統(tǒng)對導航參數(shù)誤差、傳感器測量誤差以及殘余的中、短波重力擾動分量進行估計與補償，進一步降低重力擾動對慣導解算的影響。

3重力擾動補償系統(tǒng)模型的建立

3.1狀態(tài)模型

在所設(shè)計重力擾動補償系統(tǒng)中，濾波器的狀態(tài)模型由4部分組成，包括慣性器件測量誤差模型、重力擾動狀態(tài)模型、捷聯(lián)慣導系統(tǒng)誤差模型以及陀螺/多普勒測速儀導航系統(tǒng)誤差模型。

（1）慣性器件誤差模型

對于高精度的慣性器件，可將其測量誤差建模為常值誤差與白噪聲之和，即

ε=εb+wg（16）

Δ=Δb+wa（17）

式中：ε為陀螺儀元件的測量誤差；εb和wg分別為陀螺儀測量的常值漂移和白噪聲；Δ為加速度計的測量誤差；Δb和wa分別為加速度計測量的常值偏置和白噪聲。

（2）重力擾動狀態(tài)模型

在利用EGM2008球諧模型對長波重力擾動分量進行補償后，重力擾動殘差在長波段的能量較小，因此可以利用帶阻尼的二階高斯-馬爾可夫模型來增強模型在長波段的衰減，從而與實際的重力擾動補償殘差更加貼合。該模型可表示為［15］

δg¨i（t）=－2ξω0，iδg·i（t）－ω20，iδgi（t）+qi（t）（18）

式中：i分別表示E，N，U；ξ為阻尼系數(shù)，通過調(diào)節(jié)阻尼值可以得到具有不同功率譜特性的重力擾動模型；ω0，i為馬爾科夫模型的中心頻率，可表示為2πv/di，di為重力擾動相關(guān)距離，v為載體的運動速度；qi（t）為高斯白噪聲，其方差Qi=4ξ（2πv/di）3σ2δg，i，σ2δg，i為各向重力擾動分量的方差。本文中取ξ=0.7，dE=245.4km，dN=213.3km，dU=155.5km，三向σ2δg，i由球諧模型近似計算為1069.29mGal2。

（3）捷聯(lián)慣導系統(tǒng)誤差模型

考慮重力擾動影響的慣導速度誤差方程如式（9）所示，位置誤差方程和姿態(tài)誤差方程如式（5）～（6）所示。

（4）陀螺/激光多普勒測速系統(tǒng)誤差模型

a.速度誤差方程

激光多普勒測速儀測量的是載體在其安裝坐標系下相對地面的運動速度，其解算所得載體在地理系中的速度受到姿態(tài)誤差角φD、三向安裝方位誤差角α以及標度因數(shù)誤差δk的影響，速度誤差方程為

δvnD=v^nD×φD+Cnmαv^mD+δkCnmv^mD（19）

式中：δvnD=［δvE，DδvN，DδvU，D］T為陀螺/多普勒測速儀導航系統(tǒng)解算的速度誤差。

b.位置誤差方程

對式（14）的位置更新方程等式兩邊同時求微分，可得位置誤差方程：

δL·D=δvN，DRM+h^D－V^N，D（RM+h^D）2δhDδλ·D=secL^DRN+h^DδvE，D－v^E，DsecL^D（RN+h^D）2δhD+

v^E，DsecL^DtanL^DRN+h^DδLDδh·D=δvU，D（20）

c.姿態(tài)誤差方程

對式（15）的姿態(tài)更新方程等式兩邊同時求微分，可得姿態(tài)誤差方程：

φ·nD=δωnie，D－ω^nin，D×φnD+δωnen，D－C^nb，Dδωbib（21）

選取SINS的姿態(tài)失準角φs、速度誤差δvs、位置誤差δps，陀螺/多普勒測速儀的姿態(tài)失準角φD、位置誤差δpD，陀螺的常值漂移εb，加速度計的零偏Δb，多普勒測速儀的俯仰和方位安裝誤差角α′=［αθαφ］T，標度因數(shù)誤差δk以及重力擾動及其一階導數(shù)δg和δg·作為狀態(tài)量，則狀態(tài)量共30維，即

X=［φs，δvs，δps，φD，δpD，εb，Δb，α′0，δk，δg，δg·］T（22）

將式（16）～（21）所示4部分的狀態(tài)模型聯(lián)立，可得濾波器的狀態(tài)模型：

X·（t）=F（t）X（t）+G（t）W（t）（23）

式中：F（t）為系統(tǒng)矩陣，F(xiàn)（t）=F1F206×24F3；F1與SINS誤差、陀螺/LDV導航系統(tǒng)誤差和各傳感器誤差有關(guān)，可以表示為F1=MsMD09×24；

Ms=Mφφ，sMφv，sMφp，s03×6－C^nb，s03×303×3

Mvφ，sMvv，sMvp，s03×603×3C^nb，s03×3

03×3Mpv，sMpp，s03×603×303×303×3，

Mφv，s=0－1RM+h^s01RM+h^s000tanL^sRM+h^s0，

Mφp，s=00v^N，s（RN+h^s）2

－ω^ie，ssinL^s0－v^E，s（RN+h^s）2

ω^ie，scosL^s+v^E，sRN+h^ssec2L^s0－v^E，stanL^s（RN+h^s）2，

Mvp，s=Av，s（Mω，s+Mφp，s），Mvv，s=Av，sMφv，s－Aω，s，

Mω，s=000

－ω^ie，ssinL^s00

ω^ie，scosL^s00，

Mpv，s=01RM+h^s0

secL^sRN+h^s00001，

Mpp，s=00－v^N，s（RM+h^s）2

v^E，ssecL^stanL^sRN+h^s0－v^E，ssecL^sRN+h^s000，

Mφφ，s為－ω^nin，s構(gòu)成的反對稱陣，Mvφ，s為n系下比力f^n構(gòu)成的反對稱陣，Av，s為v^ns構(gòu)成的反對稱陣，Aω，s為（2ω^nie，s+ω^nen，s）構(gòu)成的反對稱陣；

MD=03×9Mφφ，DMφp，D－C^nb，D03×3Mφk，D

03×9Mpφ，DMpp，D03×303×3Mpk，D，

Mφφ，D=Mφφ1，D－Av，DMφv，D，

Mφv，D=0－1RM+h^D0

1RN+h^D00

tanL^DRN+h^D00，

Mφφ1，D為v^N，DRM+h^D

－ω^ie，DcosL^D－v^E，DRN+h^D

－ω^ie，DsinL^D－v^E，DtanL^DRN+h^D所構(gòu)成的反對稱陣，

Mφp，D=000－ω^ie，DsinL^D00

ω^ie，DcosL^D+v^E，Dsec2L^DRN00，

Mφk，D=c23RM+hD－c21RM+hD－c22RM+hD－c13RN+hDc11RN+hDc12RN+hD

－c13tanLDRM+hDc11tanLDRM+hDc12tanLDRM+hD，

Mpφ，D=－Mpv，DAv，D，

Mpv，D=01RM+h^D0

secL^DRN+h^D00

001，

Mpp，D=000v^nE，DsecL^DtanL^DRN+h^D00000，

Mpk，D=vD－c23RM+hDc21RM+hDc22RM+hD－c13secLDRN+hDc11secLDRN+hDc12secLDRN+hD－c33c31c32，

Av，D為v^nD構(gòu)成的反對稱陣，cij為姿態(tài)矩陣C^nb，D的第i行、第j列元素；F2和F3與重力擾動有關(guān)，

F2=03×303×3－I3×303×3018×3018×3，F(xiàn)3=03×3I3×3MgMdg，

Mg=－ω20，E000－ω20，N000－ω20，U，

Mdg=－2ξω20，E000－2ξω20，N000－2ξω20，U；

W（t）為狀態(tài)噪聲，包含了陀螺和加速度計的測量白噪聲wg，wa，以及重力擾動馬爾科夫模型的驅(qū)動噪聲q，可表示為W（t）=wgwaqT；G（t）為狀態(tài)噪聲系數(shù)矩陣，G（t）=－C^nb，s03×303×3

03×3C^nb，s03×3

03×303×303×3

－C^nb，D03×303×3

015×3015×3015×3

03×303×3I3×3。

綜合以上各式可得重力擾動補償系統(tǒng)的狀態(tài)模型。

3.2量測模型

（1）位置量測模型

以SINS與陀螺/LDV系統(tǒng)解算所得位置之差Zp作為量測值，則可建立位置誤差量測模型：

Zp=p^ns－p^nD=（pn+δpns）－（pn+δpnD）=

δpns－δpnD=HpX+Vp（24）

式中：Hp=［03×6I3×303×3－I3×303×15］；Vp為虛擬位置量測噪聲。

（2）速度量測模型

以SINS與陀螺/LDV系統(tǒng)解算所得速度之差Zv作為量測值，則可建立速度誤差量測模型：

Zv=v^ns－v^nD=（vn+δvns）－（vn+δvnD）=

δvns－δvnD=HvX+Vv（25）

將式（19）陀螺/激光多普勒測速系統(tǒng)速度誤差方程代入式（25）可得速度量測矩陣的具體形式：

Hv=［03×3I3×303×3H103×9H203×6］

式中：H1為－v^nD構(gòu)成的反對稱矩陣；H2的具體形式為H2=vD－c13c11－c12－c23c21－c22－c33c31－c32；Vv為激光多普勒測速儀的速度量測噪聲。

（3）姿態(tài)量測模型

將SINS與陀螺/LDV系統(tǒng)解算所得捷聯(lián)矩陣C^nb，s轉(zhuǎn)置與C^nb，D相乘，構(gòu)造得到姿態(tài)誤差陣Cdφ：

Cdφ=C^nb，D（C^nb，s）T=（I－ΦnD）Cnb［（I－Φns）Cnb］T=

I+Φns－ΦnD（26）

式中：Φns和ΦnD分別為φns和φnD所構(gòu)成的反對稱陣。

以Cdφ的非對角線元素Zφ=φns－φnD作為量測值，可以建立姿態(tài)量測模型：

Zφ=HφX+Vφ（27）

式中：Hφ=［I3×303×6－I3×303×18］；Vφ為虛擬姿態(tài)量測噪聲。

（4）高度量測模型

將氣壓高度計測量所得hBA與慣導解算所得高度h^s之差作為量測值，則可建立高度量測模型：

Zh=HhX+Vh（28）

式中：Hh=［01×2101×27］；Vh為高度計量測噪聲。

聯(lián)立式（24）～（28）可得組合濾波器的量測模型：

Z=HX+V（29）

式中：Z為量測值；H為觀測矩陣；V為系統(tǒng)噪聲陣，具體可寫為

Z=ZpZvZφZh，H=HpHvHφHh，V=VpVvVφVh。

以上給出了重力擾動補償系統(tǒng)濾波器的狀態(tài)方程和量測方程，根據(jù)這兩組方程再加上必要的初始條件即可進行卡爾曼濾波，從而有效地估計出慣導系統(tǒng)誤差以及重力擾動矢量信息。

3.3系統(tǒng)可觀性分析

系統(tǒng)的可觀性代表由系統(tǒng)輸出或觀測得到估計狀態(tài)的可能性。系統(tǒng)能否利用量測信息有效地對狀態(tài)量進行估計，主要由系統(tǒng)的可觀性決定。因此，可觀性分析是判斷系統(tǒng)有效性的重要環(huán)節(jié)。

由于載體的運動，重力擾動補償系統(tǒng)是一個時變的系統(tǒng)，在濾波過程中狀態(tài)變量的可觀測性不恒定，給系統(tǒng)的可觀性分析帶來了困難。針對這一問題，可以利用分段定常（PWCS）的方法對系統(tǒng)進行處理［21-22］。此處以所建立的系統(tǒng)模型為基礎(chǔ)，利用分段線性定?？捎^測性分析理論對典型機動方式下系統(tǒng)的可觀性進行分析，得到系統(tǒng)可觀性矩陣的秩如表1所示，表中“#”代表不同的時段。

由表1結(jié)果可知，載體機動方式對重力擾動補償系統(tǒng)的可觀性有一定的影響。當載體以勻速直線和加速直線運動時，只有線速度或線加速度發(fā)生變化，不利于狀態(tài)量之間的解耦；當載體進行轉(zhuǎn)彎和爬升等角機動時，線速度/線加速度發(fā)生變化的同時，方向余弦和角速度也發(fā)生變化，有利于狀態(tài)量之間的解耦，因此系統(tǒng)的可觀性較好。

此外，當載體進行俯沖與爬升后，系統(tǒng)完全可觀，因此在實際應(yīng)用過程中，可先進行俯沖與爬升機動以改善系統(tǒng)的可觀性，再對重力擾動進行補償。

4仿真驗證

4.1仿真條件

運載體從30°N，120°E出發(fā)，初始速度與姿態(tài)角均為0，行駛過程中進行直線行駛、加減速、轉(zhuǎn)彎、爬坡等機動。陀螺儀的常值漂移為0.01（°）/h，量測噪聲標準差為0.005（°）/h；加速度計常值偏置為10×10-6g，零偏白噪聲為5×10-6g，三向初始位置誤差為10m，三向初始速度誤差為0.1m/s，三向初始對準誤差均為30″。LDV標度因數(shù)誤差為0.01，量測噪聲為0.01m/s；氣壓高度計量測噪聲為20m。IMU數(shù)據(jù)更新率為100Hz，LDV和氣壓高度計測量頻率均為10Hz。仿真總時間為1h，仿真步長為0.01s。

在仿真中使用的重力擾動數(shù)據(jù)由GGMplus高分地球重力數(shù)據(jù)庫提供。GGMplus是美國宇航局與德國航天局聯(lián)合研制的重力數(shù)據(jù)庫［23］，其所包含的重力場數(shù)據(jù)及其描述性統(tǒng)計結(jié)果如表2所示。

該數(shù)據(jù)庫綜合利用了重力衛(wèi)星數(shù)據(jù)、EGM重力場球諧模型與地形重力數(shù)據(jù)，能提供分辨率為7.2″的地球重力場網(wǎng)格信息，是現(xiàn)階段分辨率和精度最高的區(qū)域重力數(shù)據(jù)庫，能有效反映實際重力擾動的分布情況?；谠撃Ｐ吞峁┑闹亓_動進行仿真可有效驗證所設(shè)計方案的性能。由GGMplus數(shù)據(jù)庫得到的載體運動軌跡上重力擾動如圖7所示。

由圖7可知，仿真軌跡上的重力擾動以近似于常值的長波分量為主，并且包含局部變化較大的短波分量，與實際重力擾動特性相符。

4.2仿真結(jié)果與分析

4.2.1重力擾動的估計結(jié)果

為了驗證所提方案對重力擾動的估計效果，分別采用基于EGM球諧模型補償?shù)姆桨福▊鹘y(tǒng)方法1）、基于狀態(tài)估計補償?shù)姆桨福▊鹘y(tǒng)方法2）與所提方法對重力擾動進行估計，得到相應(yīng)垂線偏差和重力異常的估計誤差曲線，如圖8所示。

由圖8對重力擾動的估計誤差可知，直接利用EGM2008球諧模型對重力擾動的估計精度較差，對垂線偏差的估計精度約為7.17″，對重力異常的估計精度約為8.52mGal。此外，由于球諧模型僅能計算較長波段的重力擾動分量，對重力場變化的細節(jié)不敏感，因此在重力擾動變化較為劇烈時的估計精度會進一步下降。

在狀態(tài)估計補償方案中，由于重力擾動馬爾科夫模型所能覆蓋的波段有限，對長波擾動分量的描述精度較差，導致該方案對近似于常值的長波重力擾動分量估計效果不佳，即估計結(jié)果與實際值相比存在約5mGal的常值偏置，較難適用于全波段重力擾動的估計。該方案對垂線偏差的估計精度約為6.22″，對重力異常的估計精度約為7.83mGal，整體上估計精度有限。

所設(shè)計方案經(jīng)過約700s的運動后，對重力擾動的估計誤差逐漸收斂，在實際重力擾動變化較為劇烈時也能維持較好的估計效果，對垂線偏差的估計精度約為3.14″，對重力異常的估計精度約為4.20mGal。由此可知，通過利用EGM模型預先對長波重力擾動分量進行補償，可有效改善馬爾科夫模型對重力擾動殘差描述的準確度，進而提高其估計精度，從而能夠?qū)﹂L波以及中短波重力擾動分量進行有效的估計。此外，由于低階球諧模型的計算結(jié)果接近實際重力擾動的常值分量，經(jīng)過補償后可減小重力擾動殘差的常值偏置，進而減小濾波估計的初始誤差，因此當載體以同樣規(guī)律進行機動時，所設(shè)計方案對重力擾動估計誤差的收斂速度更快。

4.2.2導航解算結(jié)果

將未補償重力擾動的純慣導解算結(jié)果、利用球諧模型補償后的導航結(jié)果、利用狀態(tài)估計法補償后的導航結(jié)果以及利用所提方法補償重力擾動后的導航結(jié)果進行對比，其位置誤差曲線如圖9所示。

由圖9可知，當不對重力擾動進行補償時，導航位置誤差隨時間迅速發(fā)散。利用球諧模型對重力擾動進行補償后，能一定程度減緩導航誤差的發(fā)散。由于狀態(tài)估計法在對重力擾動估計與補償?shù)耐瑫r，能對慣導系統(tǒng)誤差進行反饋矯正，因此傳統(tǒng)方法2和所提方法的導航誤差較小。不過由于SINS和LDV航位推算均具有相同的初始誤差，而氣壓高度計只能提供高度方向的阻尼，系統(tǒng)沒有水平方向的絕對位置量測值，導致水平方向的位置誤差不收斂，東向和北向的位置誤差始終在初始值附近震蕩。

將20次仿真結(jié)果取平均，得到這幾種方法的平均導航誤差，如表3所示。

由表3結(jié)果可知，相比于單純的狀態(tài)估計法補償，所提方法通過改善重力擾動的估計與補償精度，可進一步減小重力擾動對慣導系統(tǒng)的影響，因此能夠有效提升導航精度，其水平定位精度約為3m，垂直定位精度約為0.05m，實現(xiàn)了高精度的重力擾動補償和組合導航。

5結(jié)論

實際重力擾動所覆蓋的波長范圍較寬，傳統(tǒng)方法較難同時對各波段的重力擾動分量進行高精度補償。針對這一問題，提出了一種球諧模型補償與狀態(tài)估計補償相互結(jié)合的高精度重力擾動補償方法，對長波以及中短波重力擾動分量采取兩種不同途徑進行估計。利用EGM2008球諧模型在長波段計算精度高的特點，對長波擾動分量進行精確的補償；利用馬爾科夫模型對中短波重力擾動分量描述效果較好的特點，對殘余的重力擾動分量進行高精度的狀態(tài)建模與估計，進一步提高對中短波分量的估計精度，整體上具有較好的補償效果，可實現(xiàn)高精度的重力擾動補償和組合導航。

參考文獻：

［1］WestonJL，TittertonDH.ModernInertialNavigationTechnologyandItsApplication［J］.Electronics&CommunicationEngineeringJournal，2000，12（2）：49-64.

［2］陸志東，王晶.高精度慣性導航系統(tǒng)重力補償方法［J］.航空科學技術(shù)，2016，27（8）：1-6.

LuZhidong，WangJing.GravityCompensationMethodsforHighAccuracyINS［J］.AeronauticalScience＆Technology，2016，27（8）：1-6.（inChinese）

［3］KwonJH.GravityCompensationMethodsforPrecisionINS［C］∥60thAnnualMeetingoftheInstituteofNavigation，2001：483-490.

［4］ChangLB，QinFJ，WuMP.GravityDisturbanceCompensationforInertialNavigationSystem［J］.IEEETransactionsonInstrumentationandMeasurement，2019，68（10）：3751-3765.

［5］王晶，楊功流，李湘云，等.重力擾動矢量對慣導系統(tǒng)影響誤差項指標分析［J］.中國慣性技術(shù)學報，2016，24（3）：285-290.

WangJing，YangGongliu，LiXiangyun，etal.ErrorIndicatorAnalysisforGravityDisturbingVectorsInfluenceonInertialNavigationSystem［J］.JournalofChineseInertialTechnology，2016，24（3）：285-290.（inChinese）

［6］ZhuZS，TanH，JiaY，etal.ResearchontheGravityDisturbanceCompensationTerminalforHighPrecisionPositionandOrientationSystem［J］.Sensors，2020，20（17）：4932.

［7］叢琳，趙忠，楊小步.長航時捷聯(lián)慣導重力擾動影響及補償［J］.計算機工程與應(yīng)用，2014，50（11）：251-255.

CongLin，ZhaoZhong，YangXiaobu.AnalysisandCompensationofDisturbingGravityEffectonLongEnduranceStrapDownInertialNavigationSystem［J］.ComputerEngineeringandApplications，2014，50（11）：251-255.（inChinese）

［8］WengJ，LiuJN，JiaoMX，etal.AnalysisandOnLineCompensationofGravityDisturbanceinaHighPrecisionInertialNavigationSystem［J］.GPSSolutions，2020，24（3）：1-8.

［9］PavlisNK，HolmesSA，KenyonSC，etal.TheDevelopmentandEvaluationoftheEarthGravitationalModel2008（EGM2008）［J］.JournalofGeophysicalResearch：SolidEarth，2012，117（B4）：531-535.

［10］劉雨生，樓立志.不同分辨率數(shù)字高程模型對EGM2008截斷誤差補償效果分析［J/OL］.地球物理學進展，2022.https：∥kns-cnki-net.e2.buaa.edu.cn/kcms/detail/11.2982.P.20210723.0858.012.html.

LiuYusheng，LouLizhi.AnalysisofCompensationforTruncationErrorofEGM2008byDigitalElevationModelwithDifferentResolutions［J/OL］.ProgressinGeophysics，2022.https：∥kns-cnki-net.e2.buaa.edu.cn/kcms/detail/11.2982.P.20210723.0858.012.html.（inChinese）

［11］楊金玉，張訓華，張菲菲，等.EGM2008地球重力模型數(shù)據(jù)在中國大陸地區(qū)的精度分析［J］.地球物理學進展，2012，27（4）：1298-1306.

YangJinyu，ZhangXunhua，ZhangFeifei，etal.OntheAccuracyofEGM2008EarthGravitationalModelinChineseMainland［J］.ProgressinGeophysics，2012，27（4）：1298-1306.（inChinese）

［12］WangJ，YangGL，LiXY，etal.ApplicationoftheSphericalHarmonicGravityModelinHighPrecisionInertialNavigationSystems［J］.MeasurementScienceandTechnology，2016，27（9）：095103.

［13］JekeliC.AirborneVectorGravimetryUsingPrecise，PositionAidedInertialMeasurementUnits［J］.BulletinGéodésique，1994，69（1）：1-11.

［14］AnW，XuJN，HeHY，etal.AMethodofDeflectionoftheVerticalMeasurementBasedonAttitudeDifferenceCompensation［J］.IEEESensorsJournal，2021，21（12）：13125-13136.

［15］戴東凱.基于GNSS/SINS組合的船載高精度垂線偏差測量方法研究［D］.長沙：國防科學技術(shù)大學，2016.

DaiDongkai.ResearchonHighPrecisionShipborneMeasurementMethodofDeflectionsoftheVerticalBasedonGNSS/SINSIntegration［D］.Changsha：NationalUniversityofDefenseTechnology，2016.（inChinese）

［16］趙德軍，吳曉平.空間擾動引力的譜分析［J］.海洋測繪，2005，25（2）：6-8.

ZhaoDejun，WuXiaoping.SpectralAnalysisofSpatialDisturbingGravity［J］.HydrographicSurveyingandCharting，2005，25（2）：6-8.（inChinese）

［17］JekeliC.AlgorithmsandPreliminaryExperienceswiththeLN93andLN100forAirborneVectorGravimetry［R］.DefenseTechnicalInformationCenter，1998.

［18］MoritzH.GeodeticReferenceSystem1980［J］.BulletinGéodésique，1980，54（3）：395-405.

［19］鐵俊波，吳美平，蔡劭琨，等.基于EGM2008重力場球諧模型的水平重力擾動計算方法［J］.中國慣性技術(shù)學報，2017，25（5）：624-629.

TieJunbo，WuMeiping，CaiShaokun，etal.GravityDisturbanceCalculationMethodBasedonEarthGravitationalModel2008［J］.JournalofChineseInertialTechnology，2017，25（5）：624-629.（inChinese）

［20］WangQ，GaoCF，ZhouJ，etal.TwoDimensionalLaserDopplerVelocimeterandItsIntegratedNavigationwithaStrapdownInertialNavigationSystem［J］.AppliedOptics，2018，57（13）：3334-3339.

［21］周衛(wèi)東，蔡佳楠，孫龍，等.GPS/SINS超緊組合導航系統(tǒng)可觀測性分析［J］.北京航空航天大學學報，2013，39（9）：1157-1162.

ZhouWeidong，CaiJianan，SunLong，etal.ObservabilityAnalysisofGPS/SINSUltraTightlyCoupledSystem［J］.JournalofBeijingUniversityofAeronauticsandAstronautics，2013，39（9）：1157-1162.（inChinese）

［22］肖佳敏，朱鋒，張小紅.GNSS/SINS松組合系統(tǒng)的可觀測性分析［J］.導航定位學報，2018，6（4）：35-41.

XiaoJiamin，ZhuFeng，ZhangXiaohong.ObservabilityAnalysisonLooselyCoupledGNSS/SINSIntegrationSystem［J］.JournalofNavigationandPositioning，2018，6（4）：35-41.（inChinese）

［23］HirtC，ClaessensS，F(xiàn)echerT，etal.NewUltrahighResolutionPictureofEarthsGravityField［J］.GeophysicalResearchLetters，2013，40（16）：4279-4283.

[HJ*3][HJ][JZ（]AnAccurateGravityDisturbanceCompensationMethod

BasedonSphericalHarmonicModelandMultiSensorFusion

LiuYuxin1，WangXinlong1*，WangXun2，GaoWenning3，HuXiaodong4

（1.SchoolofAstronautics，BeihangUniversity，Beijing100083，China;

2.BeijingInstituteofAutomaticControlEquipment，Beijing100074，China;

3.BeijingInstituteofSatelliteInformationEngineering，Beijing100095，China；

4.AVICXianFlightAutomaticControlResearchInstitute，Xian710065，China）

Abstract：Withtheimprovementoftheaccuracyrequirementsofhighprecisionlongenduranceinertialnavigationsystem，thegravitydisturbancehasbecomethemainerrorsourceofinertialnavigationsystem.Isthekeyfactortoimprovethenavigationaccuracywhetheritcanbeeffectivelycompensated.Thetraditionalcompensationmethodbasedonsphericalharmonicmodelcannotreflectthedetailedinformationoftheearthsgravityfield，andthushaspoorcompensationresultformediumandshortwavegravitydisturbancecomponents.Inthetraditionalstateestimationmethod，theMarkovstatemodelhaspooraccuracyindescribingthelongwavegravitydisturbancecomponent，sothecompensationresultofthelongwavecomponentispoor.Theabovemethodscannotcompensatetheactualgravitydisturbancethathaswidebandwithhighprecision.Tosolvethisproblem，ahighprecisiongravitydisturbancecompensationmethodbasedonsphericalharmonicmodelandmultisensorinformationfusionisdesignedinthispaper.Ontheonehand，themethodcompensatesthelongwavegravitydisturbancecomponentbyusingthehighcalculationaccuracyofthelowordersphericalharmonicmodelinthelongwaveband.Ontheotherhand，usingstrapdowninertialnavigation/laserDopplervelocimetry/barometricaltimetertoformafullyautonomousintegratednavigationsystem，theresidualmediumandshortwavegravitydisturbancecomponentisestablishedasahighprecisionMarkovmodel，soastorealizethestateestimationandcompensationofmediumandshortwavecomponent.Thesimulationresultsshowthattheproposedgravitydisturbancecompensationmethodcaneffectivelyimprovetheestimationeffectofgravitydisturbance，realizehighprecisiongravitydisturbancecompensationandhavehighnavigationaccuracy.

Keywords：strapdowninertialnavigationsystem;gravitydisturbancecompensation;sphericalharmonicgravitymodel;stateestimation;integratednavigation

收稿日期：2022-03-01

基金項目：國家自然科學基金項目（61673040）；重點基礎(chǔ)研究項目（2020-JCJQ-ZD-136-12）；航空科學基金項目（20170151002）；天地一體化信息技術(shù)國家重點實驗室基金項目（2015-SGIIT-KFJJ-DH-01）

作者簡介：劉宇鑫（1999-），男，北京人，碩士研究生。

*通信作者：王新龍（1969-），男，陜西渭南人，教授。