[發明專利]一種慣導系統誤差仿真及修理輔助分析方法

在審

• 申請號: 202110705955.7
• 申請日: 2021-06-24
• 公開（公告）號: CN113432622A
• 公開（公告）日: 2021-09-24
• 發明（設計）人: 覃剛;范云鵬;姚宇;林曉彬;董海靜;董金珠
• 申請（專利權）人: 中國船舶重工集團公司第七0七研究所
• 主分類號: G01C25/00
• 分類號: G01C25/00
• 代理公司: 天津盛理知識產權代理有限公司 12209
• 代理人: 王雨晴

1.一種慣導系統誤差仿真及修理輔助分析方法，其特征在于：包括以下步驟：

步驟1、建立慣性導航系統誤差方程；

步驟2、根據步驟1所建立的慣性導航系統誤差方程，設置初始條件，對慣導系統誤差進行仿真及修理輔助分析，并根據仿真結果對影響慣導系統精度的決定性因素進行分析。

2.根據權利要求1所述的一種慣導系統誤差仿真及修理輔助分析方法，其特征在于：所述步驟1的具體步驟包括：

(1)建立靜基座誤差方程：

(2)寫成狀態方程矩陣形式并進行拉氏變換得到：

3.根據權利要求1所述的一種慣導系統誤差仿真及修理輔助分析方法，其特征在于：所述步驟2的具體步驟包括：

(1)設置初始條件

初始條件：設陀螺的常值漂移為0.1°/h，加速度計的常值零偏誤差為0.0001g，當地緯度39°；起始條件誤差為：速度誤差0.1m/s，位置誤差為0.0005°，水平姿態誤差為20角秒，方位姿態誤差為5角分。運行時間t＝24小時；

(2)對慣導系統誤差進行仿真及修理輔助分析，推導出陀螺漂移引起的系統誤差，并畫出誤差傳播曲線；推導出加速度計零偏引起的系統誤差，并畫出誤差傳播曲線；推導出初始條件誤差引起的系統誤差，并畫出誤差傳播曲線；

(3)根據仿真結果對影響慣導系統精度的決定性因素進行分析。

4.根據權利要求3所述的一種慣導系統誤差仿真及修理輔助分析方法，其特征在于：所述步驟2第(2)步的具體步驟包括：

①建立系數陣為：

用Matlab求特征矩陣C＝(sI-F)^-1有：

C₁₁＝(R*s)/(R*s^2+g)

C₁₂＝0

C₁₃＝(R*g*s*wie*sin(L))/((s^2+wie^2)*(R*s^2+g))

C₁₄＝(R*g*s*wie*sin(L))/((s^2+wie^2)*(R*s^2+g))

C₁₅＝-(R*g*(s^2-wie^2*sin(L)^2+wie^2))/((s^2+wie^2)*(R*s^2+g))

C₁₆＝-(R*g*wie^2*sin(2*L))/(2*(s^2+wie^2)*(R*s^2+g))

C₂₁＝0

C₂₂＝(R*s)/(R*s^2+g)

C₂₃＝-(R*g*wie^2)/((s^2+wie^2)*(R*s^2+g))

C₂₄＝(R*g*s^2)/((s^2+wie^2)*(R*s^2+g))

C₂₅＝(R*g*s*wie*sin(L))/((s^2+wie^2)*(R*s^2+g))

C₂₆＝-(R*g*s*wie*cos(L))/((s^2+wie^2)*(R*s^2+g))

C₃₁＝0

C₃₂＝1/(R*s^2+g)

C₃₃＝(R*s^3+(R*wie^2+g)*s)/((s^2+wie^2)*(R*s^2+g))

C₃₄＝(g*s)/((s^2+wie^2)*(R*s^2+g))

C₃₅＝(g*wie*sin(L))/((s^2+wie^2)*(R*s^2+g))

C₃₆＝-(g*wie*cos(L))/((s^2+wie^2)*(R*s^2+g))

C₄₁＝0

C₄₂＝-1/(R*s^2+g)

C₄₃＝-(R*s*wie^2)/((s^2+wie^2)*(R*s^2+g))

C₄₄＝(R*s^3)/((s^2+wie^2)*(R*s^2+g))

C₄₅＝(R*s^2*wie*sin(L))/((s^2+wie^2)*(R*s^2+g))

C₄₆＝-(R*s^2*wie*cos(L))/((s^2+wie^2)*(R*s^2+g))

C₅₁＝1/(R*s^2+g)

C₅₂＝0

C₅₃＝-(R*s^2*wie*sin(L))/((s^2+wie^2)*(R*s^2+g))

C₅₄＝-(R*s^2*wie*sin(L))/((s^2+wie^2)*(R*s^2+g))

C₅₅＝(R*s*(s^2-wie^2*sin(L)^2+wie^2))/((s^2+wie^2)*(R*s^2+g))

C₅₆＝(R*s*wie^2*sin(2*L))/(2*(s^2+wie^2)*(R*s^2+g))

C₆₁＝tan(L)/(R*s^2+g)

C₆₂＝0

C₆₃＝(wie*(R*s^2*cos(L)^2+g))/(cos(L)*(s^2+wie^2)*(R*s^2+g))

C₆₄＝(wie*(R*s^2*cos(L)^2+g))/(cos(L)*(s^2+wie^2)*(R*s^2+g))

C₆₅＝-(s*sin(L)*(g-R*wie^2*cos(L)^2))/(cos(L)*(s^2+wie^2)*(R*s^2+g))

C₆₆＝(R*s^3+(R*wie^2*sin(L)^2+g)*s)/((s^2+wie^2)*(R*s^2+g))

各式中，常數值g表示重力加速度，ws表示舒勒角頻率，wie表示地球自轉角速度，R表示地球半徑，L表示當地緯度；

Dx表示東向加速度計常值零偏，Dy表示北向加速度計常值零偏，Ex表示東向陀螺漂移，Ey表示北向陀螺漂移，Ez表示方位陀螺漂移；

Dx0表示初始東向速度誤差，Dy0表示初始北向速度誤差，L0表示初始緯度位置誤差，N0表示初始經度位置誤差，Qx0表示初始東向姿態誤差，Qy0表示初始北向姿態誤差，Qz0表示初始方位姿態誤差；

dVx表示系統東向速度誤差，dVy表示系統北向速度誤差；dL表示緯度誤差，dN表示經度誤差；Qx、Qy、Qz表示平臺誤差；

②計算東向陀螺漂移εx引起的系統誤差及誤差傳播曲線；

dVxs＝(Ex*R*g*wie*sin(L))/((s^2+wie^2)*(R*s^2+g))

dVxt＝(Ex*R*g*sin(t*wie)*sin(L))/(g-R*wie^2)-(Ex*R^(3/2)*g^(1/2)*wie*sin((g^(1/2)*t)/R^(1/2))*sin(L))/(g-R*wie^2)

dVys＝(Ex*R*g*s)/((s^2+wie^2)*(R*s^2+g))

dVyt＝(Ex*R*g*cos(t*wie))/(g-R*wie^2)-(Ex*R*g*cos((g^(1/2)*t)/R^(1/2)))/(g-R*wie^2)

dLs＝(Ex*g)/((s^2+wie^2)*(R*s^2+g))

dLt＝(Ex*g*sin(t*wie))/(wie*(g-R*wie^2))-(Ex*R^(1/2)*g^(1/2)*sin((g^(1/2)*t)/R^(1/2)))/(g-R*wie^2)

dNs＝(Ex*g*wie*sin(L))/(cos(L)*(s^2+wie^2)*(R*s^2+g))

dNt＝((Ex*g*sin(t*wie)*sin(L))/(g-R*wie^2)-(Ex*R^(1/2)*g^(1/2)*wie*sin((g^(1/2)*t)/R^(1/2))*sin(L))/(g-R*wie^2))/cos(L)

Qxs＝(Ex*R*s^2)/((s^2+wie^2)*(R*s^2+g))

Qxt＝(Ex*R^(1/2)*g^(1/2)*sin((g^(1/2)*t)/R^(1/2)))/(g-R*wie^2)-(Ex*R*wie*sin(t*wie))/(g-R*wie^2)

Qys＝-(Ex*R*s*wie*sin(L))/((s^2+wie^2)*(R*s^2+g))

Qyt＝(Ex*R*wie*cos((g^(1/2)*t)/R^(1/2))*sin(L))/(g-R*wie^2)-(Ex*R*wie*cos(t*wie)*sin(L))/(g-R*wie^2)

Qzs＝(Ex*wie*(R*s^2*cos(L)^2+g))/(s*cos(L)*(s^2+wie^2)*(R*s^2+g))

Qzt＝-((Ex*cos(t*wie)*(g-R*wie^2*cos(L)^2))/(wie*(g-R*wie^2))-Ex/wie+(Ex*R*wie*cos((g^(1/2)*t)/R^(1/2))*(cos(L)-1)*(cos(L)+1))/(g-R*wie^2))/cos(L)

③計算北向陀螺漂移εy引起的系統誤差及誤差傳播曲線；

dVxs＝-(Ey*R*g*(s^2-wie^2*sin(L)^2+wie^2))/(s*(s^2+wie^2)*(R*s^2+g))

dVxt＝Ey*R*(sin(L)-1)*(sin(L)+1)+(Ey*R*cos((g^(1/2)*t)/R^(1/2))*(R*wie^2*sin(L)^2-R*wie^2+g))/(g-R*wie^2)-(Ey*R*g*cos(t*wie)*sin(L)^2)/(g-R*wie^2)

dVys＝(Ey*R*g*wie*sin(L))/((s^2+wie^2)*(R*s^2+g))

dVyt＝(Ey*R*g*sin(t*wie)*sin(L))/(g-R*wie^2)-(Ey*R^(3/2)*g^(1/2)*wie*sin((g^(1/2)*t)/R^(1/2))*sin(L))/(g-R*wie^2)

dLs＝(Ey*g*wie*sin(L))/(s*(s^2+wie^2)*(R*s^2+g))

dLt＝(Ey*sin(L))/wie+(Ey*R*wie*cos((g^(1/2)*t)/R^(1/2))*sin(L))/(g-R*wie^2)-(Ey*g*cos(t*wie)*sin(L))/(wie*(g-R*wie^2))

dNs＝-(Ey*g*(s^2-wie^2*sin(L)^2+wie^2))/(s*cos(L)*(s^2+wie^2)*(R*s^2+g))

dNt＝(Ey*(sin(L)-1)*(sin(L)+1)+(Ey*cos((g^(1/2)*t)/R^(1/2))*(R*wie^2*sin(L)^2-R*wie^2+g))/(g-R*wie^2)-(Ey*g*cos(t*wie)*sin(L)^2)/(g-R*wie^2))/cos(L)

Qxs＝(Ey*R*s*wie*sin(L))/((s^2+wie^2)*(R*s^2+g))

Qxt＝(Ey*R*wie*cos(t*wie)*sin(L))/(g-R*wie^2)-(Ey*R*wie*cos((g^(1/2)*t)/R^(1/2))*sin(L))/(g-R*wie^2)

Qys＝(Ey*R*(s^2-wie^2*sin(L)^2+wie^2))/((s^2+wie^2)*(R*s^2+g))

Qyt＝(Ey*R^(1/2)*sin((g^(1/2)*t)/R^(1/2))*(R*wie^2*sin(L)^2-R*wie^2+g))/(g^(1/2)*(g-R*wie^2))-(Ey*R*wie*sin(t*wie)*sin(L)^2)/(g-R*wie^2)

Qzs＝-(Ey*sin(L)*(g-R*wie^2*cos(L)^2))/(cos(L)*(s^2+wie^2)*(R*s^2+g))

Qzt＝-((Ey*sin(t*wie)*sin(L)*(g-R*wie^2*cos(L)^2))/(wie*(g-R*wie^2))-(Ey*R^(1/2)*sin((g^(1/2)*t)/R^(1/2))*sin(L)*(g-R*wie^2*cos(L)^2))/(g^(1/2)*(g-R*wie^2)))/cos(L)

④計算方位陀螺漂移εz引起的系統誤差及誤差傳播曲線；

dVxs＝-(Ez*R*g*wie^2*sin(2*L))/(s*(R*s^2+g)*(2*s^2+2*wie^2))

dVxt＝(Ez*R*g*sin(2*L)*cos(t*wie))/(2*g-2*R*wie^2)-(Ez*R^2*wie^2*sin(2*L)*cos((g^(1/2)*t)/R^(1/2)))/(2*g-2*R*wie^2)-(Ez*R*sin(2*L))/2

dVys＝-(Ez*R*g*wie*cos(L))/((s^2+wie^2)*(R*s^2+g))

dVyt＝(Ez*R^(3/2)*g^(1/2)*wie*cos(L)*sin((g^(1/2)*t)/R^(1/2)))/(g-R*wie^2)-(Ez*R*g*cos(L)*sin(t*wie))/(g-R*wie^2)

dLs＝-(Ez*g*wie*cos(L))/(s*(s^2+wie^2)*(R*s^2+g))

dLt＝(Ez*g*cos(t*wie)*cos(L))/(wie*(g-R*wie^2))-(Ez*R*wie*cos((g^(1/2)*t)/R^(1/2))*cos(L))/(g-R*wie^2)-(Ez*cos(L))/wie

dNs＝-(Ez*g*wie^2*sin(2*L))/(s*cos(L)*(R*s^2+g)*(2*s^2+2*wie^2))

dNt＝-(Ez*sin(2*L)-(Ez*g*sin(2*L)*cos(t*wie))/(g-R*wie^2)+(Ez*R*wie^2*sin(2*L)*cos((g^(1/2)*t)/R^(1/2)))/(g-R*wie^2))/(2*cos(L))

Qxs＝-(Ez*R*s*wie*cos(L))/((s^2+wie^2)*(R*s^2+g))

Qxt＝(Ez*R*wie*cos((g^(1/2)*t)/R^(1/2))*cos(L))/(g-R*wie^2)-(Ez*R*wie*cos(t*wie)*cos(L))/(g-R*wie^2)

Qys＝(Ez*R*wie^2*sin(2*L))/((R*s^2+g)*(2*s^2+2*wie^2))

Qyt＝(Ez*R*wie*sin(2*L)*sin(t*wie))/(2*g-2*R*wie^2)-(Ez*R^(3/2)*wie^2*sin(2*L)*sin((g^(1/2)*t)/R^(1/2)))/(g^(1/2)*(2*g-2*R*wie^2))

Qzs＝(Ez*(R*s^3+(R*wie^2*sin(L)^2+g)*s))/(s*(s^2+wie^2)*(R*s^2+g))

Qzt＝(Ez*sin(t*wie)*(R*wie^2*sin(L)^2-R*wie^2+g))/(wie*(g-R*wie^2))-(Ez*R^(3/2)*wie^2*sin((g^(1/2)*t)/R^(1/2))*sin(L)^2)/(g^(1/2)*(g-R*wie^2))

⑤計算東向加速度計零偏引起的系統誤差及誤差傳播曲線；

dVxs＝(Dx*R)/(R*s^2+g)

dVxt＝(Dx*R^(1/2)*sin((g^(1/2)*t)/R^(1/2)))/g^(1/2)

dVys＝0

dVyt＝0

dLs＝0

dLt＝0

dNs＝Dx/(cos(L)*(R*s^2+g))

dNt＝(Dx*sin((g^(1/2)*t)/R^(1/2)))/(R^(1/2)*g^(1/2)*cos(L))

Qxs＝0

Qxt＝0

Qys＝Dx/(s*(R*s^2+g))

Qyt＝Dx/g-(Dx*cos((g^(1/2)*t)/R^(1/2)))/g

Qzs＝(Dx*tan(L))/(s*(R*s^2+g))

Qzt＝(Dx*tan(L))/g-(Dx*cos((g^(1/2)*t)/R^(1/2))*tan(L))/g

⑥計算初始方位姿態誤差Φ_z0引起的系統誤差及誤差傳播曲線；

dVxs＝-(Qz0*R*g*wie^2*sin(2*L))/((R*s^2+g)*(2*s^2+2*wie^2))

dVxt＝(Qz0*R^(3/2)*g^(1/2)*wie^2*sin(2*L)*sin((g^(1/2)*t)/R^(1/2)))/(2*g-2*R*wie^2)-(Qz0*R*g*wie*sin(2*L)*sin(t*wie))/(2*g-2*R*wie^2)

dVys＝-(Qz0*R*g*s*wie*cos(L))/((s^2+wie^2)*(R*s^2+g))

dVyt＝(Qz0*R*g*wie*cos((g^(1/2)*t)/R^(1/2))*cos(L))/(g-R*wie^2)-(Qz0*R*g*wie*cos(t*wie)*cos(L))/(g-R*wie^2)

dLs＝-(Qz0*g*wie*cos(L))/((s^2+wie^2)*(R*s^2+g))

dLt＝(Qz0*R^(1/2)*g^(1/2)*wie*cos(L)*sin((g^(1/2)*t)/R^(1/2)))/(g-R*wie^2)-(Qz0*g*cos(L)*sin(t*wie))/(g-R*wie^2)

dNs＝-(Qz0*g*wie^2*sin(2*L))/(cos(L)*(R*s^2+g)*(2*s^2+2*wie^2))

dNt＝-((Qz0*g*wie*sin(2*L)*sin(t*wie))/(g-R*wie^2)-(Qz0*R^(1/2)*g^(1/2)*wie^2*sin(2*L)*sin((g^(1/2)*t)/R^(1/2)))/(g-R*wie^2))/(2*cos(L))

Qxs＝-(Qz0*R*s^2*wie*cos(L))/((s^2+wie^2)*(R*s^2+g))

Qxt＝(Qz0*R*wie^2*cos(L)*sin(t*wie))/(g-R*wie^2)-(Qz0*R^(1/2)*g^(1/2)*wie*cos(L)*sin((g^(1/2)*t)/R^(1/2)))/(g-R*wie^2)

Qys＝(Qz0*R*s*wie^2*sin(2*L))/((R*s^2+g)*(2*s^2+2*wie^2))

Qyt＝(Qz0*R*wie^2*sin(2*L)*cos(t*wie))/(2*g-2*R*wie^2)-(Qz0*R*wie^2*sin(2*L)*cos((g^(1/2)*t)/R^(1/2)))/(2*g-2*R*wie^2)

Qzs＝(Qz0*(R*s^3+(R*wie^2*sin(L)^2+g)*s))/((s^2+wie^2)*(R*s^2+g))

Qzt＝(Qz0*cos(t*wie)*(R*wie^2*sin(L)^2-R*wie^2+g))/(g-R*wie^2)-(Qz0*R*wie^2*cos((g^(1/2)*t)/R^(1/2))*sin(L)^2)/(g-R*wie^2)。