[發明專利]ECEF坐標系下2-D雷達誤差配準兩步算法

無效

• 申請號: 201110124564.2
• 申請日: 2011-05-16
• 公開（公告）號: CN102305927A
• 公開（公告）日: 2012-01-04
• 發明（設計）人: 王國宏;陳壘;劉德浩;賈舒宜
• 申請（專利權）人: 中國人民解放軍海軍航空工程學院
• 主分類號: G01S7/40
• 分類號: G01S7/40
• 代理公司: 暫無信息
• 代理人: 暫無信息

1.一種基于ECEF坐標系的兩坐標雷達誤差配準方法，其特征在于包括以下步驟：

步驟1：雷達i(i取A、B，下同)將k時刻本地地理坐標：緯度Ls_i，經度Rs_i，高度H_Si和對目標的量測值：距離r_i，方位角θ_i上報到融合中心；

步驟2：融合中心使用雷達A、B上報的信息進行處理，得到等效量測方程；

在融合中心數據處理計算機中執行以下步驟：

(1)分別計算雷達i在ECEF坐標系中的坐標和旋轉矩陣：

xis(k)=(Ci(k)+Hsi(k))cos(Lsi(k))cos(λsi(k))yis(k)=(Ci(k)+Hsi(k))cos(Lsi(k))sin(λsi(k))zis(k)=(Ci(k)(1-e2)+Hsi(k))sin(Lsi(k))---(1)]]>

其中，Eq＝6378137.0，b＝6356752.3142；

Ti(k)=-sin(Rsi(k))-sin(Lsi(k))cos(Rsi(k))cos(Lsi(k))cos(Rsi(k))cos(Rsi(k))-sin(Lsi(k))sin(Rsi(k))cos(Lsi(k))sin(Rsi(k))0cos(Lsi(k))sin(Lsi(k))---(2)]]>

Di1(k)=Ti(k)-sin(θi(k))-ri(k)cos(θi(k))-cos(θi(k))ri(k)sin(θi(k))00---(3)]]>

Z1(k)=xBS(k)-xAS(k)yBS(k)-yAS(k)0+TB(k)rB(k)sin(θB(k))rB(k)cos(θB(k))0-TA(k)rA(k)sin(θA(k))rA(k)cos(θA(k))0---(4)]]>

F₁(k)＝G₁(k)＝[D_A1(k)-D_B1(k)]???????????????????????(5)

步驟3：Kalman濾波初始化

(1)Kalman濾波值初始化

β^(k|k)=(0,0,0,0)T---(6)]]>

為系統偏差估計值，為1×4向量，分別表示Δr_A(k)、Δθ_A(k)、Δr_B(k)、Δθ_B(k)，這里k＝1，即在第1時刻進行初始化，上標“^T”表示矩陣轉置；

(2)Kalman濾波協方差矩陣P(k|k)初始化

P(k|k)=10700001000010700001---(7)]]>

步驟4：Kalman濾波更新

(1)狀態預測值為：

β^(k+1|k)=β^(k|k)---(8)]]>

(2)量測預測值為：

Z^1(k+1|k)=F1(k+1)β^(k+1|k)---(9)]]>

(3)新息協方差矩陣為：

S(k+1)=F1(k+1)P(k|k)F1T(k+1)+G1(k+1)RG1T(k+1)---(10)]]>

其中，R表示雷達量測噪聲方差陣，為4×4對角陣，分別為雷達A距離量測噪聲方差，雷達A方位角量測噪聲方差，雷達B距離量測噪聲方差，雷達B方位角量測噪聲方差；

(4)增益矩陣為：

K(k+1)=P(k|k)F1T(k+1)S-1(k+1)---(11)]]>

(5)狀態估計更新

β^(k+1|k+1)=β^(k+1|k)+K(k+1)(Z1(k+1)-Z^1(k+1|k))---(12)]]>

(6)狀態估計協方差更新

P(k+1|k+1)=(I(4)-K(k+1)F1(k+1))P(k|k)(I(4)-K(k+1)F1(k+1))T]]>

+K(k+1)G1(k+1)RG1T(k+1)KT(k+1)---(13)]]>

其中，I(4)為4×4單位陣；

步驟5：對計算結果進行修正

(1)設置高度值；

H＝90000m

RnA(k+1)=-12H2rA(k+1)-18H4rA3(k+1)---(14)]]>

RnB(k+1)=-12H2rB(k+1)-18H4rB3(k+1)]]>

(2)生成等效方程

Di2(k+1)=Ti0-Rni(k+1)cos(θi(k+1))0Rni(k+1)sin(θi(k+1))00---(15)]]>

Z2(k+1)=00zBs-zAs+TB(k+1)RnB(k+1)sin(θB(k+1))RnB(k+1)cos(θB(k+1))H-TA(k+1)RnA(k+1)sin(θA(k+1))RnA(k+1)cos(θA(k+1))H---(16)]]>

F₂(k+1)＝[D_A2(k+1)-D_B2(k+1)]??????????????????(17)

(3)計算修正量

Z(k+1)=Z2(k+1)-F2(k+1)×β^(k+1)---(18)]]>

a1=F1(1,1)+F2(1,1)F1(1,3)+F2(1,3)F1(2,1)+F2(2,1)F1(2,3)+F2(2,3)---(19)]]>

b1=Z(1,1)F1(1,3)+F2(1,3)Z(2,1)F1(2,3)+F2(2,3)---(20)]]>

c1=F1(1,1)+F2(1,1)Z(1,1)F1(2,1)+F2(2,1)Z(2,1)---(21)]]>

其中，F_i(u，v)表示矩陣F_i(k+1)的第u行第v列元素；

Δβ^(1)=det(b1)/det(a1)---(22)]]>

Δβ^(3)=det(c1)/det(a1)---(23)]]>

其中，det(A)表示求矩陣A的行列式；

(4)對解的可信度進行判斷

如果c，則否則

如果則否則

其中，表示第i行元素；

步驟6：令k＝k+1，重復執行步驟1、2、4、5進行下一周期的系統偏差估計，直至雷達關機，其中，步驟4中的狀態估計初始值和狀態估計協方差為上一次Kalman濾波得到的狀態估計值和狀態估計協方差值，而非步驟5得到的修正后的狀態估計值和狀態估計協方差值。