機載SAR慣性導航誤差快速估計方法*

王冠勇，王風飛，張磊，胡慶榮

(1.北京無線電測量研究所，北京 100854；2.西安電子科技大學 雷達信號處理國家重點實驗室，陜西 西安 710071)

0 引言

合成孔徑雷達(synthetic aperture radar, SAR)[1]是全天時對地觀測的有效手段，廣泛運用于星載、機載、彈載等平臺，其高分辨成像為環境監測、災害評估、軍事打擊等應用提供技術支撐[2]。相比于對SAR圖像質量的評估，SAR定位精度[3]也是衡量SAR系統性能的一項重要指標。對于無人機載SAR[4-5]，受成本和載荷的限制，無法配備精度較高的機載慣導系統，因此慣導誤差成為引起圖像定位誤差主要因素之一，能否準確而又快速估計慣導誤差是實現無人機載SAR高定位精度實時成像的關鍵因素。

傳統的機載SAR定位算法有多項式對地定位法、共線方程對地定位法[6]和距離-多普勒(range-Doppler, R-D)模型對地定位法[7]等，其中多項式對地定位法僅對圖像作變形處理，需要大量地面控制點(ground control point, GCP)，共線方程對地定位法和R-D模型對地定位法通過反演載機真實位置進行定位修正，已被廣泛應用[8-11]。但此類方法僅局限于SAR成像結果的后處理，無法實現在飛行中對慣導誤差的實時估計。

針對這個問題，本文研究了一種結合子圖像后向投影(back projection, BP)成像的慣導誤差快速估計方法。通過分析慣導誤差在BP成像模型下引起定位誤差，可以建立定位誤差觀測方程。并且由于GCP的精確地理坐標是已知的，因此可以對以各GCP為中心的子圖像網格進行BP成像，從而高效獲取GCP的定位誤差。與傳統的R-D模型對地定位法相比，本文算法可以快速解算慣導運動誤差，并且在相同均值和方差的觀測噪聲條件下具有更小的克拉美羅界(Cramer-Rao bound, CRB)。下面將對算法原理進行詳細說明。

1 BP成像模型定位誤差分析

1.1 成像模型

BP算法成像幾何模型如圖1所示。成像幾何坐標系中x軸為沿速度方向即OA方向，y軸垂直于速度方向即BO′方向，z軸指向地面即OB方向。點P為待成像目標點，載機與點P存在相對運動關系，(xp,yp,zp)為P點在成像幾何坐標系下以O點為坐標原點在x,y,z軸的坐標。成像平面為地理坐標網格，可以將圖像的定位誤差在成像平面內向x軸和y軸方向投影，分別表示為Δx和Δy。為了簡化分析過程，假設慣導在一個合成孔徑時間內存在恒定的位置誤差為ΔP=(ΔX,ΔY,ΔZ)T，以及恒定的速度誤差為Δv=(Δvx,Δvy,Δvz)T，那么可以通過等效斜距歷史顯式求解目標點在BP成像平面中的定位誤差。下面將分別分析各種慣導誤差在BP成像模型下的定位誤差。

圖1 SAR成像幾何模型Fig 1 Geometric model of SAR imaging

1.2 x軸方向位置誤差

對于場景內任意目標點P，理想情況下雷達天線相位中心到P點的瞬時斜距歷史Rp表示為

(1)

式中：X=vtm為載機沿航向運動軌跡；v為載機理想速度；tm為方位慢時間。在慣導存在x軸方向位置誤差ΔX的情況下，真實的斜距歷史表達式為

(2)

觀察式(2)可知，該式等效于理想天線相位中心到點目標P′的斜距歷史，其中P′點與P點在X軸方向相距Δx。因此可以判斷慣導存在x軸方向位置誤差ΔX的情況下，成像結果在x軸方向發生平移，x軸方向定位誤差Δx為

Δx=ΔX

.

(3)

1.3 y軸方向位置誤差

在慣導存在y軸方向位置誤差ΔY的情況下，真實的斜距歷史表達式為

(4)

式(4)等效于理想天線相位中心到點目標P′的斜距歷史，點P′與P點相比在y軸方向相距ΔY，則成像結果在y軸方向定位誤差為

Δy=ΔY.

(5)

1.4 z軸方向位置誤差

當慣導存在z軸方向位置誤差ΔZ時，真實的斜距歷史表達式為

(6)

為了表示z軸方向定位誤差在成像網格坐標系內的投影，可以將式(6)變換為

(7)

式中：Δy表示z軸方向誤差ΔZ在圖像網格中沿y軸方向的等效誤差，其表達式為

(8)

將式(8)在ΔZ=0附近進行Taylor展開，保留一次項系數，得到成像結果y軸方向的定位誤差為

(9)

1.5 x軸方向速度誤差

當慣導存在x軸方向速度誤差Δvx時，天線相位中心到目標點P的斜距歷史為

(10)

為了用等效點目標的形式表示式(10)，可將式(10)近似等效為

(11)

式中：Δx，Δy為等效目標點在x軸方向和y軸方向的定位誤差，表達式為

(12)

(13)

將式(13)在Δvz=0附近進行Taylor展開，保留一次項，得到y軸方向定位誤差近似表達式為

(14)

下面來驗證式(11)相對式(10)的等效精度，將式(10)在tm=0附近進行Taylor展開，保留二階多項式，得

(15)

式中：

(16)

同樣地，將式(11)在tm=0附近進行Taylor展開，保留二階多項式，得

(17)

式中：

(18)

1.6 y軸方向速度誤差

當慣導存在y軸方向速度誤差Δvy時，天線相位中心到目標點P的斜距歷史為

(19)

為了顯式表示運動誤差下目標點的定位誤差，可將式(19)近似等效為

(20)

式中：Δx，Δy為等效目標點的在x軸方向和y軸方向的定位誤差，表達式為

(21)

(22)

將式(22)在Δvy=0附近進行Taylor展開，保留一次項，得到y軸方向定位誤差近似表達式為

(23)

1.7 z軸方向速度誤差

當慣導存在z軸方向速度誤差Δvz時，天線相位中心到目標點P的斜距歷史為

(24)

為了顯式表示運動誤差下目標點的定位誤差，可將式(24)近似等效為

(25)

式中:Δx,Δy為等效目標點的在x軸方向和y軸方向定位誤差，表達式為

(26)

(27)

將式(27)在Δvz=0附近進行Taylor展開，保留一次項系數，得到y軸方向定位誤差近似表達式為

(28)

綜上所述，慣導位置和速度誤差會使目標點在BP成像模型中的沿航向(x軸方向)和垂直航向(y軸方向)產生定位誤差，二者關系表達式如表1所示，因此可以利用這個關系通過觀測一定數量已知GCP的定位誤差來估計慣導位置和速度誤差。

表1 BP模型定位誤差與慣導位置和速度誤差的關系

2 慣導誤差估計算法

2.1 算法原理

根據上述分析可知，在BP成像算法中，慣導的位置和速度誤差會導致成像結果中目標點在x軸方向和y軸方向分別產生定位誤差。因此通過測量足夠數量的GCP在SAR圖像中的定位誤差，就可以反向估計慣導的位置和速度誤差。由于目標點定位誤差與慣導誤差近似成線性關系，因此可以利用最小二乘法[12]進行快速估計。

此外，需要說明的是，基于BP成像模型的估計方式相比傳統R-D模型估計方法還具有一個明顯的優勢。傳統R-D法是一種SAR圖像后處理方法，需要將整幅圖像成像完畢后，在圖像中尋找GCP進行慣導誤差估計和修正，這種方法并不能運用于對實時性要求較高的SAR成像系統中，例如彈載SAR、無人機載SAR等。而子圖像BP成像模型可以將已知GCP地理坐標作為中心點建立子圖像成像網格進行快速成像，避免了對其余絕大部分區域的成像運算，縮短了GCP定位誤差的獲取時間。需要說明的是，通過結合現有的快速BP算法[13-15]和并行運算方式[16]，可以進一步加速成像運算。

假設共有m個GCP，各子圖像成像結果中目標點定位誤差觀測方程如式(29)所示，該線性方程建立了慣導誤差與圖像定位誤差之間的關系。

G=HL+n

，

(29)

式中：

G=(Δx1,Δy1,…,Δxm,Δym)T

,

(30)

L=(ΔX,ΔY,ΔZ,Δvx,Δvy,Δvz)T

,

(31)

(32)

E(nnT)=Cn

.

(33)

在式(29)中，G為觀測的目標點定位誤差向量。假設對于下角標K，ΔxK和ΔrK分別表示第K個 GCP在x方向和y方向上的定位誤差。P為待估計的慣導誤差向量，ΔX,ΔY,ΔZ,Δvx,Δvy,Δvz分別表示慣導位置和速度誤差在成像幾何坐標系x軸、y軸和z軸方向上的投影。n為觀測噪聲向量，自相關矩陣表示為Cn。H∈2m×6為觀測矩陣，xK,yK,zK分別表示第K個GCP在成像幾何坐標系中的坐標。根據最小二乘估計法，誤差向量P的估計值為

(34)

2.2 估計性能分析

(35)

(36)

(37)

對比式(36)和式(37)可知，在Gauss-Markov定理條件下，即誤差向量的各獨立分布的高斯隨機變量均具有零均值和相同方差σ2(Cn=σ2I，I是單位矩陣)時，最小二乘估計的均方誤差矩陣等于CRB。

3 仿真實驗

實驗在給定的一組位置和速度誤差下進行，分別用基于傳統距離多普勒成像的R-D模型估計方法與基于BP成像模型的估計方法進行對比。在不同標準差的觀測噪聲條件下分別進行10 000次Monte Carlo實驗，比較6個待估計參數在2種估計算法下的CRB和估計均方誤差隨觀測噪聲的變化曲線。仿真參數如表2所示。

表2 仿真實驗參數

假設BP圖像中每個點目標x軸方向與y軸方向的定位誤差觀測相互獨立，R-D圖像中斜距和多普勒的觀測相互獨立，觀測噪聲均滿足均值為0的高斯分布。在BP成像模型的估計方法中，觀測噪聲協方差矩陣為

(38)

式中：σg為目標點位置(或斜距)觀測的標準差；矩陣I2M是大小為2M×2M的單位矩陣，M為地面控制點的數量。對于距離—多普勒估計方法，觀測噪聲協方差矩陣為

(39)

式中：σf為目標點多普勒頻率觀測標準差，多普勒頻率觀測精度受斜距觀測精度影響，σf與σg近似滿足如下關系式：

(40)

仿真實驗以觀測噪聲標準差σg為變量，比較不同σg下2種算法估計的均方誤差。σg的變化范圍為0～3 m，根據式(40)的對應關系，σf的變化范圍為0～2.54 Hz。

設待估計向量為P1=(3,5,8,0.6,0.5,-0.5)T，第1組地面控制點坐標如表3所示，通過實驗得到6個待估參量在2種估計方法下的均方誤差與CRB如圖2所示。圖2a)～f)分別為在傳統R-D模型估計法和BP成像模型的估計法下，通過10 000次Monte Carlo實驗得到對比結果，其中包含慣導位置和速度誤差分別在x,y和z方向估計的均方誤差與CRB理論值的對應結果。圖2結果說明，在相同觀測噪聲均方誤差條件下，BP成像模型的慣導誤差估計方法與傳統R-D模型慣導誤差估計方法相比可以得到更小的估計均方誤差。

表3 地面控制點坐標

圖2 2種估計方法下慣導位置、速度誤差估計均方誤差與相應CRB的對比Fig. 2 Comparison of mean square error results of inertial navigation position and velocity error estimations with corresponding CRBs under two estimation methods

4 結束語

針對機載SAR低精度慣導誤差估計問題，本文提出了一種基于子圖像BP定位的快速估計方法。通過BP成像方式，可以利用并行運算實現各GCP所在子圖像的快速成像，加速慣導誤差估計過程。此外，在等量觀測噪聲條件下，本文算法與傳統R-D模型估計法相比具有更小的估計均方誤差，仿真實驗進行了驗證。

[1] CUMMING I,WONG F.Digital Processing of Synthetic Aperture Radar Data:Algorithm and Implementation[M].Norwood,MA:Artech House,2005.

[2] 邢濤，胡慶榮，李軍,等.毫米波高分辨SAR成像算法性能分析[J].現代防御技術,2015,43(1):81-86. XING Tao,HU Qing-rong,LI Jun,et al.Analysis of Millimeter Wave High Resolution SAR Imaging Algorithm[J].Modern Defence Technology,2015,43(1):81-86.

[3] CURLANDER J C.Location of Pixels in Space-Borne SAR Imagery[J].IEEE Transactions on Geoscience and Remote Sensing,1982,20(3):359-364.

[4] ZHANG Lei,XING Meng-dao,QIAO Zhi-jun.Wavenumber-Domain Autofocusing for Highly Squinted UAV SAR Imagery[J].IEEE Sensor Journal,2012,12(5):1574-1588.

[5] ZHANG Lei,QIAO Zhi-jun,XING Meng-dao,et al.A Robust Motion Compensation Approach for UAV SAR Imagery[J]. IEEE Transactions on Geoscience and Remote Sensing,2012,50(8):3202-3218.

[6] KONECNY G,SCHUHR W.Reliability of Radar Image Data[C]∥Proceedings of 16th ISPRS Congress,Kyoto,1988:92-101.

[7] LEBERA F.Radargammetric Imaging Processing[M].MA:Artech House,1990.

[8] TANNOUS I,PIKEROEB B.Parametric Modeling of Spaceborne SAR Image Geometry Application:SEASAT/SPOT Image Registration[J].Photogrammetric Engineering and Remote Sensing,1994,60(6):755-766.

[9] DAVID S,FRANCESCO H,ERICH M,et al.Geometric and Radiometric Calibration of RadarSat Images[C]∥Proceedings of Geomatics in the Era of RadarSat,Ottawa,1997:24-30.

[10] JOHNEN H,LAUKNES L,GUNERIUSSEN T.Geocoding of Fast-Delivery ERS-1 SAR Image Mode Product Using DEM Data[J].International Journal of Remote Sensing,1995,16(11):1957-1968.

[11] VASSILLIA K,CHRISTOS L.SAR Geocoding Method for Evaluation Geodetic Coordinates and Improving Indirect Geocoding Accuracy[C]∥International Symposium on Remote Sensing,2002:268-278.

[12] GROEN P D.An Introduction to Total Least Squares[J].Nieuw Archief voor Wiskunde,1996,14(2):237-254.

[13] YEGULALP A F.Fast Backprojection Algorithm for Synthetic Aperture Radar[C]∥International Proceedings of Radar Conference,Waltham,MA,USA,Apr.20-22,1999：60-65.

[14] ULANDER L M H,HELLSTEN H,STENSTROM G..Synthetic Aperture Radar Processing Using Fast Factorized Back-Projection[J].IEEE Transactions on Aerospace Electronics System,2003,39(3):760-776.

[15] ZHANG Lei,LI Hao-lin,QIAO Zhi-jun,et al.A Fast BP Algorithm with Wavenumber Spectrum Fusion for High-Resolution Spotlight SAR Imaging[J].IEEE Geoscience and Remote Sensing Letters,2014,11(9):1460-1464.

[16] CHRISTOPHE E,MICHEL J,INFLADA J.Remote Sensing Processing:From Multicore to GPU[J].IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing,2011,4(3):643-652.