王冠勇，王風(fēng)飛，張磊，胡慶榮

(1.北京無線電測(cè)量研究所，北京 100854；2.西安電子科技大學(xué) 雷達(dá)信號(hào)處理國家重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710071)

0 引言

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

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

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

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

1.1 成像模型

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

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

1.2 x軸方向位置誤差

對(duì)于場(chǎng)景內(nèi)任意目標(biāo)點(diǎn)P，理想情況下雷達(dá)天線相位中心到P點(diǎn)的瞬時(shí)斜距歷史Rp表示為

(1)

式中：X=vtm為載機(jī)沿航向運(yùn)動(dòng)軌跡；v為載機(jī)理想速度；tm為方位慢時(shí)間。在慣導(dǎo)存在x軸方向位置誤差ΔX的情況下，真實(shí)的斜距歷史表達(dá)式為

(2)

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

Δx=ΔX

.

(3)

1.3 y軸方向位置誤差

在慣導(dǎo)存在y軸方向位置誤差ΔY的情況下，真實(shí)的斜距歷史表達(dá)式為

(4)

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

Δy=ΔY.

(5)

1.4 z軸方向位置誤差

當(dāng)慣導(dǎo)存在z軸方向位置誤差ΔZ時(shí)，真實(shí)的斜距歷史表達(dá)式為

(6)

為了表示z軸方向定位誤差在成像網(wǎng)格坐標(biāo)系內(nèi)的投影，可以將式(6)變換為

(7)

式中：Δy表示z軸方向誤差ΔZ在圖像網(wǎng)格中沿y軸方向的等效誤差，其表達(dá)式為

(8)

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

(9)

1.5 x軸方向速度誤差

當(dāng)慣導(dǎo)存在x軸方向速度誤差Δvx時(shí)，天線相位中心到目標(biāo)點(diǎn)P的斜距歷史為

(10)

為了用等效點(diǎn)目標(biāo)的形式表示式(10)，可將式(10)近似等效為

(11)

式中：Δx，Δy為等效目標(biāo)點(diǎn)在x軸方向和y軸方向的定位誤差，表達(dá)式為

(12)

(13)

將式(13)在Δvz=0附近進(jìn)行Taylor展開，保留一次項(xiàng)，得到y(tǒng)軸方向定位誤差近似表達(dá)式為

(14)

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

(15)

式中：

(16)

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

(17)

式中：

(18)

1.6 y軸方向速度誤差

當(dāng)慣導(dǎo)存在y軸方向速度誤差Δvy時(shí)，天線相位中心到目標(biāo)點(diǎn)P的斜距歷史為

(19)

為了顯式表示運(yùn)動(dòng)誤差下目標(biāo)點(diǎn)的定位誤差，可將式(19)近似等效為

(20)

式中：Δx，Δy為等效目標(biāo)點(diǎn)的在x軸方向和y軸方向的定位誤差，表達(dá)式為

(21)

(22)

將式(22)在Δvy=0附近進(jìn)行Taylor展開，保留一次項(xiàng)，得到y(tǒng)軸方向定位誤差近似表達(dá)式為

(23)

1.7 z軸方向速度誤差

當(dāng)慣導(dǎo)存在z軸方向速度誤差Δvz時(shí)，天線相位中心到目標(biāo)點(diǎn)P的斜距歷史為

(24)

為了顯式表示運(yùn)動(dòng)誤差下目標(biāo)點(diǎn)的定位誤差，可將式(24)近似等效為

(25)

式中:Δx,Δy為等效目標(biāo)點(diǎn)的在x軸方向和y軸方向定位誤差，表達(dá)式為

(26)

(27)

將式(27)在Δvz=0附近進(jìn)行Taylor展開，保留一次項(xiàng)系數(shù)，得到y(tǒng)軸方向定位誤差近似表達(dá)式為

(28)

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

表1 BP模型定位誤差與慣導(dǎo)位置和速度誤差的關(guān)系

2 慣導(dǎo)誤差估計(jì)算法

2.1 算法原理

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

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

假設(shè)共有m個(gè)GCP，各子圖像成像結(jié)果中目標(biāo)點(diǎn)定位誤差觀測(cè)方程如式(29)所示，該線性方程建立了慣導(dǎo)誤差與圖像定位誤差之間的關(guān)系。

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為觀測(cè)的目標(biāo)點(diǎn)定位誤差向量。假設(shè)對(duì)于下角標(biāo)K，ΔxK和ΔrK分別表示第K個(gè) GCP在x方向和y方向上的定位誤差。P為待估計(jì)的慣導(dǎo)誤差向量，ΔX,ΔY,ΔZ,Δvx,Δvy,Δvz分別表示慣導(dǎo)位置和速度誤差在成像幾何坐標(biāo)系x軸、y軸和z軸方向上的投影。n為觀測(cè)噪聲向量，自相關(guān)矩陣表示為Cn。H∈2m×6為觀測(cè)矩陣，xK,yK,zK分別表示第K個(gè)GCP在成像幾何坐標(biāo)系中的坐標(biāo)。根據(jù)最小二乘估計(jì)法，誤差向量P的估計(jì)值為

(34)

2.2 估計(jì)性能分析

(35)

(36)

(37)

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

3 仿真實(shí)驗(yàn)

實(shí)驗(yàn)在給定的一組位置和速度誤差下進(jìn)行，分別用基于傳統(tǒng)距離多普勒成像的R-D模型估計(jì)方法與基于BP成像模型的估計(jì)方法進(jìn)行對(duì)比。在不同標(biāo)準(zhǔn)差的觀測(cè)噪聲條件下分別進(jìn)行10 000次Monte Carlo實(shí)驗(yàn)，比較6個(gè)待估計(jì)參數(shù)在2種估計(jì)算法下的CRB和估計(jì)均方誤差隨觀測(cè)噪聲的變化曲線。仿真參數(shù)如表2所示。

表2 仿真實(shí)驗(yàn)參數(shù)

假設(shè)BP圖像中每個(gè)點(diǎn)目標(biāo)x軸方向與y軸方向的定位誤差觀測(cè)相互獨(dú)立，R-D圖像中斜距和多普勒的觀測(cè)相互獨(dú)立，觀測(cè)噪聲均滿足均值為0的高斯分布。在BP成像模型的估計(jì)方法中，觀測(cè)噪聲協(xié)方差矩陣為

(38)

式中：σg為目標(biāo)點(diǎn)位置(或斜距)觀測(cè)的標(biāo)準(zhǔn)差；矩陣I2M是大小為2M×2M的單位矩陣，M為地面控制點(diǎn)的數(shù)量。對(duì)于距離—多普勒估計(jì)方法，觀測(cè)噪聲協(xié)方差矩陣為

(39)

式中：σf為目標(biāo)點(diǎn)多普勒頻率觀測(cè)標(biāo)準(zhǔn)差，多普勒頻率觀測(cè)精度受斜距觀測(cè)精度影響，σf與σg近似滿足如下關(guān)系式：

(40)

仿真實(shí)驗(yàn)以觀測(cè)噪聲標(biāo)準(zhǔn)差σg為變量，比較不同σg下2種算法估計(jì)的均方誤差。σg的變化范圍為0～3 m，根據(jù)式(40)的對(duì)應(yīng)關(guān)系，σf的變化范圍為0～2.54 Hz。

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

表3 地面控制點(diǎn)坐標(biāo)

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

4 結(jié)束語

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

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

[2] 邢濤，胡慶榮，李軍,等.毫米波高分辨SAR成像算法性能分析[J].現(xiàn)代防御技術(shù),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.