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

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

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

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

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

關鍵詞：捷聯慣導系統；重力擾動補償；重力場球諧模型；狀態估計；組合導航

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

文獻標識碼：A

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

DOI：10.12132/ISSN.1673-5048.2022.0038

0引言

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

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

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

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

針對以上問題，本文設計了一種球諧模型補償與狀態估計補償相互結合的重力擾動補償方法，在缺乏區域重力網格數據庫時能完全自主地實現重力擾動的估計與補償，并對該方法的性能進行了驗證。

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

1.1重力擾動成因

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

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

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

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

δg=g－γ0（2）

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

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

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

由此可見，重力擾動會引起實際重力與理論重力在當地地理坐標系下三個方向的誤差。

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

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

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

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

δ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為地球自轉角速度；δωnie，s為地球自轉角速度的計算誤差；ω^nen，s為導航系相對地球系的轉動角速度；δωnie，s為導航系相對地球系的轉動角速度的計算誤差；δvns=［δvE，sδvN，sδvU，s］為速度誤差；δpns=［δLsδλsδhs］T為位置誤差；φs為姿態失準角；δg0為理論重力加速度的計算誤差；δfb為比力測量誤差；δωbib為陀螺儀的測量誤差。

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

gn=gn0+δgn0（7）

僅采用理論重力加速度公式對比力測量值進行補償，會導致計算結果中存在一定的重力擾動殘差，進而引起額外的加速度解算誤差，此時加速度解算誤差δ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為加速度計測量誤差。

受重力擾動影響的加速度誤差經過兩次積分后分別引起速度誤差δ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）引起額外的姿態解算誤差。姿態誤差又進一步通過式（4）～（5）造成位置誤差和速度誤差，形成閉環回路，持續對慣導系統的位置、速度和姿態解算產生影響。重力擾動對捷聯慣導系統的影響機理與加速度計零偏類似，其所造成的慣導解算誤差隨時間積累。由以上分析可知，重力擾動是影響高精度長航時慣導系統性能的關鍵因素之一。

2球諧模型與狀態估計相結合的補償方案設計

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

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

2.1球諧模型補償

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

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

λres=πaN（12）

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

2.2狀態估計補償

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

2.2.1陀螺/激光多普勒測速系統

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

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

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

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

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。

與捷聯慣導系統類似，利用陀螺儀測量所得角速度ω～bib以及由位置和速度計算所得的ω^nie，D和ω^nen，D，可對姿態矩陣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構成的反對稱矩陣。

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

2.2.2狀態估計補償系統

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

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

2.3重力擾動補償方案設計

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

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

3重力擾動補償系統模型的建立

3.1狀態模型

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

（1）慣性器件誤差模型

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

ε=εb+wg（16）

Δ=Δb+wa（17）

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

（2）重力擾動狀態模型

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

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

式中：i分別表示E，N，U；ξ為阻尼系數，通過調節阻尼值可以得到具有不同功率譜特性的重力擾動模型；ω0，i為馬爾科夫模型的中心頻率，可表示為2πv/di，di為重力擾動相關距離，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）捷聯慣導系統誤差模型

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

（4）陀螺/激光多普勒測速系統誤差模型

a.速度誤差方程

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

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

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

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.姿態誤差方程

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

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

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

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

將式（16）～（21）所示4部分的狀態模型聯立，可得濾波器的狀態模型：

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

式中：F（t）為系統矩陣，F（t）=F1F206×24F3；F1與SINS誤差、陀螺/LDV導航系統誤差和各傳感器誤差有關，可以表示為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構成的反對稱陣，Mvφ，s為n系下比力f^n構成的反對稱陣，Av，s為v^ns構成的反對稱陣，Aω，s為（2ω^nie，s+ω^nen，s）構成的反對稱陣；

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所構成的反對稱陣，

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構成的反對稱陣，cij為姿態矩陣C^nb，D的第i行、第j列元素；F2和F3與重力擾動有關，

F2=03×303×3－I3×303×3018×3018×3，F3=03×3I3×3MgMdg，

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

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

W（t）為狀態噪聲，包含了陀螺和加速度計的測量白噪聲wg，wa，以及重力擾動馬爾科夫模型的驅動噪聲q，可表示為W（t）=wgwaqT；G（t）為狀態噪聲系數矩陣，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。

綜合以上各式可得重力擾動補償系統的狀態模型。

3.2量測模型

（1）位置量測模型

以SINS與陀螺/LDV系統解算所得位置之差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系統解算所得速度之差Zv作為量測值，則可建立速度誤差量測模型：

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

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

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

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

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

（3）姿態量測模型

將SINS與陀螺/LDV系統解算所得捷聯矩陣C^nb，s轉置與C^nb，D相乘，構造得到姿態誤差陣Cdφ：

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

I+Φns－ΦnD（26）

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

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

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

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

（4）高度量測模型

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

Zh=HhX+Vh（28）

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

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

Z=HX+V（29）

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

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

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

3.3系統可觀性分析

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

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

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

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

4仿真驗證

4.1仿真條件

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

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

該數據庫綜合利用了重力衛星數據、EGM重力場球諧模型與地形重力數據，能提供分辨率為7.2″的地球重力場網格信息，是現階段分辨率和精度最高的區域重力數據庫，能有效反映實際重力擾動的分布情況。基于該模型提供的重力擾動進行仿真可有效驗證所設計方案的性能。由GGMplus數據庫得到的載體運動軌跡上重力擾動如圖7所示。

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

4.2仿真結果與分析

4.2.1重力擾動的估計結果

為了驗證所提方案對重力擾動的估計效果，分別采用基于EGM球諧模型補償的方案（傳統方法1）、基于狀態估計補償的方案（傳統方法2）與所提方法對重力擾動進行估計，得到相應垂線偏差和重力異常的估計誤差曲線，如圖8所示。

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

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

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

4.2.2導航解算結果

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

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

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

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

5結論

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

參考文獻：

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

［2］陸志東，王晶.高精度慣性導航系統重力補償方法［J］.航空科學技術，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］王晶，楊功流，李湘云，等.重力擾動矢量對慣導系統影響誤差項指標分析［J］.中國慣性技術學報，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］叢琳，趙忠，楊小步.長航時捷聯慣導重力擾動影響及補償［J］.計算機工程與應用，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］劉雨生，樓立志.不同分辨率數字高程模型對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地球重力模型數據在中國大陸地區的精度分析［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］.長沙：國防科學技術大學，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］.中國慣性技術學報，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］周衛東，蔡佳楠，孫龍，等.GPS/SINS超緊組合導航系統可觀測性分析［J］.北京航空航天大學學報，2013，39（9）：1157-1162.

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

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

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

［23］HirtC，ClaessensS，FecherT，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）；重點基礎研究項目（2020-JCJQ-ZD-136-12）；航空科學基金項目（20170151002）；天地一體化信息技術國家重點實驗室基金項目（2015-SGIIT-KFJJ-DH-01）

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

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