[發明專利]基于INS和GPS組合的垂線偏差動態測量裝置及方法

有效

• 申請號: 201410305314.2
• 申請日: 2014-06-30
• 公開（公告）號: CN104061945B
• 公開（公告）日: 2016-11-30
• 發明（設計）人: 王省書;戴東凱;戰德軍;秦石喬;黃宗升;鄭佳興;吳偉;胡春生
• 申請（專利權）人: 中國人民解放軍國防科學技術大學
• 主分類號: G01C25/00
• 分類號: G01C25/00;G06F19/00
• 代理公司: 國防科技大學專利服務中心 43202
• 代理人: 李振

1.基于INS和GPS組合的垂線偏差動態測量裝置，其特征在于：該裝置由慣性導航系統(1)、GPS天線(2)、GPS接收機(3)和數據處理計算機(4)組成，所述慣性導航系統(1)、GPS天線(2)以及GPS接收機(3)組成INS/GPS姿態測量子系統，所述慣性導航系統(1)包含了三個正交安裝的激光陀螺，稱為激光陀螺組合體，所述激光陀螺組合體、GPS天線(2)以及GPS接收機(3)組成LGU/GPS姿態測量子系統，三個激光陀螺均與GPS接收機(3)通信；所述慣性導航系統(1)、GPS天線(2)、GPS接收機(3)固聯安裝于測量載體(5)上，所述GPS天線(2)與所述GPS接收機(3)通信，所述慣性導航系統(1)、GPS接收機(3)通過數據線與數據處理計算機(4)連接，慣性導航系統(1)和GPS接收機(3)的測量數據通過數據線傳輸到數據處理計算機(4)中，在數據處理計算機(4)中完成垂線偏差的解算。

2.一種如權利要求1所述INS和GPS組合的垂線偏差動態測量裝置，其特征在于：所述測量載體(5)為測量船或測量車。

3.一種如權利要求1所述INS和GPS組合的垂線偏差動態測量裝置，其特征在于：所述INS采用單軸旋轉式激光陀螺慣性導航系統。

4.一種利用如權利要求1所述測量裝置動態測量垂線偏差的方法，其特征在于該方法包括以下步驟：

1)構建INS/GPS姿態測量子系統，并進行8小時以上對準：

啟動INS/GPS姿態測量子系統，進行8小時以上對準，在對準過程中INS/GPS姿態測量子系統實現組合姿態測量，在整個測量過程中INS/GPS姿態測量子系統連續輸出INS坐標系相對于計算導航坐標系的姿態矩陣對準過程中無垂線偏差數據輸出；

2)INS/GPS對準結束后，啟動LGU/GPS姿態測量子系統，利用第一步中INS/GPS對準后輸出的姿態矩陣對LGU/GPS進行姿態矩陣的初始化：

記t_i為第i個測量采樣點對應的時刻，1≤i≤N，N為整個測量過程的采樣點總數，t_i時刻INS/GPS輸出的姿態矩陣為LGU/GPS輸出的姿態為上標n表示真實導航坐標系，啟動LGU/GPS姿態測量子系統的時刻記為t₀，LGU/GPS姿態矩陣的初始化方法為令Cb(0)n(0)=Cb(0)n′(0);]]>

3)LGU/GPS姿態測量子系統進行姿態更新：

LGU/GPS姿態測量子系統進行姿態更新，同時數據處理計算機存儲全部測量采樣時刻INS/GPS輸出的姿態矩陣和LGU/GPS輸出的姿態矩陣

4)利用全部測量過程中t_i時刻INS/GPS輸出的姿態矩陣和LGU/GPS輸出的姿態矩陣計算坐標旋轉矩陣并計算其輸出的相應的姿態角之差：

的計算方法為：

Cn′(i)n(i)=Cb(i)n(i)·[Cb(i)n′(i)]T---(1)]]>

其中[]^T表示矩陣的轉置，符號·表示矩陣乘法；

將任意t時刻計算得到的n′系到n系的姿態矩陣簡記為任意t時刻INS/GPS與LGU/GPS輸出的三個姿態角之差為：ΔΦ_E，ΔΦ_N，ΔΦ_U，其中，下標E、N、U分別表示東向、北向、天向分量，ΔΦ_E、ΔΦ_N和ΔΦ_U由下式計算得到：

ΔΦE=Cn′n(2,3)ΔΦN=Cn′n(3,1)ΔΦU=Cn′n(1,2)---(2)]]>

其中表示矩陣的第2行，第3列的元素，表示矩陣的第3行，第1列的元素，表示矩陣的第1行，第2列的元素；

5)以INS/GPS與LGU/GPS輸出的三個姿態角之差ΔΦ_E，ΔΦ_N，ΔΦ_U為觀測量，通過建立垂線偏差測量的觀測方程和狀態方程，在全球重力模型的輔助下，利用Kalman濾波的方法提取垂線偏差：

記t時刻測量載體所在位置真實的東向垂線偏差為η，北向垂線偏差為ξ，由全球重力模型計算得到的東向和北向垂線偏差值分別記為由全球重力模型計算得到的東向和北向垂線偏差誤差分別記為δη和δξ，并存在下式的關系：

η=η^+δηξ=ξ^+δξ---(3)]]>

記INS/GPS輸出的相對于n′系的姿態誤差角為v_E，v_N，v_U，LGU/GPS輸出的相對于n系的姿態誤差角為φ_E、φ_N、φ_U；具體實施方法如下：

5.1)建立垂線偏差測量的狀態方程：

選取垂線偏差測量系統的狀態變量為φ_E、φ_N、φ_U、δη、δξ、ε_U，其中ε_U為激光陀螺的等效天向零偏，分別對上述狀態變量進行動態建模，φ_E、φ_N、φ_U滿足如下微分方程：

φ·E=-ωiecosL·φU+ωiesinL·φNφ·N=-ωiesinL·φEφ·U=ωiecosL·φE-ϵU---(4)]]>

其中，ω_ie為地球自轉角速度，L為測量點的地理緯度；

ε_U建模為隨機常值模型，則有：

ϵ·U=0---(5)]]>

全球重力模型的東向和北向垂線偏差誤差δη、δξ的統計模型分別由以下兩式給出：

x·Eδη·=01-ω02-2ζω0xEδη+0wE---(6)]]>

x·Nδξ·=01-ω02-2ζω0xNδξ+0wN---(7)]]>

其中x_E、x_N為中間狀態變量，ω₀為該統計模型的固有頻率，ω₀與測量載體的運動速度V之間存在固定關系ω₀＝2π×V/10000，ζ為該統計模型的阻尼系數，w_E和w_N分別為東向和北向垂線偏差統計模型的過程噪聲，w_E～(0,q_E)，w_N～(0,q_N)，即w_E和w_N分別服從均值為0的高斯分布，方差分別為q_E和q_N；

選擇垂線偏差測量系統的狀態空間矢量為：

x＝[φ_E,φ_N,φ_U,ε_U,x_E,x_N,δη,δξ]^T(8)

將式(4)-(7)寫為統一的狀態方程的形式為：

x·=F·x+w---(9)]]>

其中，矩陣F=0ωiesinL-ωiecosL00000-ωiesinL0000000ωiecosL00-100000000000000000010000000010000-ω020-2ζω0000000-ω020-2ζω0,]]>

w＝[0,0,0,0,0,0,w_E,w_N]^T，式(9)即為垂線偏差測量系統的狀態方程，過程噪聲w～(0,Q)，即w服從均值為0，方差為Q的高斯分布，Q=O6×6O6×2O2×6diag[qE,qN],]]>其中O_6×6、O_2×6、O_6×2分別表示6行6列、2行6列、6行2列的零矩陣，diag[q_E,q_N]表示以q_E和q_N為對角元素的對角矩陣；

將狀態方程(9)離散化得到：

x_i＝M_i/i-1·x_i-1+w_i(10)

x_i表示x在t_i時刻的采樣值,M_i/i-1為t_i-1到t_i時刻的狀態轉移矩陣，M_i/i-1≈I_8×8+Δt·F_i，其中I_8×8為8維的單位矩陣，F_i為矩陣F在t_i時刻的采樣值，Δt為采樣間隔，w_i服從均值為0，方差為Q_i的高斯分布，Q_i＝Δt²·Q；

5.2)建立垂線偏差測量的觀測方程：

以INS/GPS和LGU/GPS輸出的姿態角之差ΔΦ_E，ΔΦ_N為觀測量，建立如下觀測方程：

ΔΦE=ξ-φE+vEΔΦN=-η-φN+vN---(11)]]>

將式(3)代入式(11)，整理后得到：

ΔΦE-ξ^=δξ-φE+vEΔΦN+η^=-δη-φN+vN---(12)]]>

選取新的觀測矢量為：

y=[ΔΦE-ξ^,ΔΦN+η^]T---(13)]]>

將式(12)重新寫為：

y＝H·x+v(14)

其中矩陣H=-100000010-10000-10,]]>v＝[v_E,v_N]^T；

式(14)即為垂線偏差測量系統的觀測方程，INS/GPS的姿態誤差角v_E、v_N分別為東向和北向觀測分量的觀測噪聲，v_E～(0,r_E)，v_N～(0,r_N)，即v_E和v_N分別服從均值為0的高斯分布，方差分別為r_E和r_N，觀測噪聲v～(0,R)，即v服從均值為0，方差為R的高斯分布，觀測噪聲方差矩陣r=rE00rN;]]>

將觀測方程離散化得到：

y_i＝H_i·x_i+v_i(15)

y_i、v_i、H_i分別為y、v、H在t_i時刻的采樣值，v_i服從均值為0，方差為R_i的高斯分布，R_i＝Δt²·R，H_i＝H；

5.3)根據式(10)和(15)所描述的垂線偏差測量系統離散形式的狀態方程和觀測方程，利用Kalman濾波算法對狀態矢量進行估計：

Kalman濾波算法的迭代算法如下：

x^i/i-1=Mi/i-1·x^i-1Pi/i-1=Mi/i-1·Pi-1Mi/i-1T+QiKi=Pi/i-1·HiT·[Hi·Pi/i-1·HiT+Ri]-1Pi=(I-Ki·Hi)·Pi/i-1·[I-Ki·Hi]-1+Ki·Pi/i-1·KiTx^i=x^i/i-1+Ki·(yi-Hi·x^i/i-1)---(16)]]>

其中為x_i的估計值，為x_i的一步預測值，P_i為x_i的估計協方差矩陣，P_i/i-1為P_i的一步預測值，K_i為濾波增益，[]^-1符號表示矩陣的求逆運算；

t_i時刻東向和北向垂線偏差誤差δη、δξ的估計值分別為由t_i時刻的狀態矢量估計值給出：

δη^i=x^i(7)δξ^i=x^i(8)---(17)]]>

其中，分別表示狀態矢量估計值的第7和第8個元素；

5.4)計算垂線偏差測量值：

t_i時刻東向和北向垂線偏差測量值η_i、ξ_i分別由式(18)、(19)計算得到：

ηi=η^i+δη^i---(18)]]>

ξi=ξ^i+δξ^i---(19)]]>

其中分別為全球重力模型計算得到的t_i時刻測量點上東向和北向垂線偏差的值。

5.一種利用如權利要求4所述測量裝置動態測量垂線偏差的方法，其特征在于：所述全球重力模型采用的是EGM2008全球重力模型。