劉 輝，李葛爽，徐 青，樓良盛

1. 華北水利水電大學(xué)資源與環(huán)境學(xué)院，河南 鄭州 450046； 2. 鄭州市交通規(guī)劃勘察設(shè)計(jì)研究院，河南 鄭州 450001； 3. 信息工程大學(xué)，河南 鄭州 450001； 4. 西安測(cè)繪研究所，陜西 西安 710054

1996年，文獻(xiàn)[1]首次提出了通過(guò)跨航向線性陣列天線實(shí)現(xiàn)三維成像的概念。這個(gè)里程碑式的想法標(biāo)志著陣列SAR技術(shù)的誕生，拓展了SAR系統(tǒng)下視工作能力，有效解決了傳統(tǒng)SAR體制機(jī)底盲區(qū)、幾何失真和左右模糊等問(wèn)題[2-5]。該技術(shù)在跨航向上稀疏布置收發(fā)分置的陣列天線[6-8]，采用天底觀測(cè)的方式，分別通過(guò)合成孔徑、脈沖壓縮和波束形成技術(shù)實(shí)現(xiàn)觀測(cè)對(duì)象方位向、高程向和跨航向的三維分辨。目前陣列SAR技術(shù)的研究主要集中在基于平臺(tái)勻速直線運(yùn)動(dòng)假設(shè)的成像算法上[9-15]。但實(shí)際上，平臺(tái)運(yùn)動(dòng)中不可避免會(huì)產(chǎn)生位置偏移和姿態(tài)變化，從而造成天線相位中心的偏移，影響高分辨率陣列SAR三維成像質(zhì)量。因此為了保證成像所需的相位精度，分析平臺(tái)位移、姿態(tài)變化的誤差補(bǔ)償方法具有重要的意義。

在國(guó)外，法國(guó)ONERA實(shí)驗(yàn)室研制了基于BUSARD滑翔機(jī)的DRIVE Project，并利用POS系統(tǒng)獲取的位置和姿態(tài)信息輔助完成側(cè)視和正視成像試驗(yàn)，但沒(méi)有公布其運(yùn)動(dòng)補(bǔ)償方案和結(jié)果[16-17]。德國(guó)FGAN-FHR的ARTINO系統(tǒng)分別利用DGPS、IMU組合設(shè)備和激光/CCD裝置，分析姿態(tài)角誤差和機(jī)翼振動(dòng)對(duì)成像的影響，但前者也未公布補(bǔ)償方案和結(jié)果，后者雖提出了補(bǔ)償方法但未考慮目標(biāo)空變性影響[18-21]。2010年后國(guó)外學(xué)者再無(wú)公開發(fā)表過(guò)陣列SAR相關(guān)研究成果，目前公開發(fā)表的文獻(xiàn)主要來(lái)自國(guó)內(nèi)學(xué)者。2009年，文獻(xiàn)[22]提出了下視聚束三維SAR PGA補(bǔ)償方法，但相位誤差非空變假設(shè)不符合下視陣列SAR觀測(cè)幾何；2012年，文獻(xiàn)[23]研究了機(jī)載三維SAR多通道聯(lián)合自聚焦補(bǔ)償方法，但忽略了目標(biāo)空變性影響；文獻(xiàn)[24]分析了相位中心偏差對(duì)機(jī)載陣列天線下視3D-SAR成像的影響。2015年，文獻(xiàn)[25—26]研究了機(jī)載三維SAR平動(dòng)誤差、偏航角誤差補(bǔ)償方法。但這些研究都是從信號(hào)角度分析的，與測(cè)繪關(guān)聯(lián)甚少。

本文從測(cè)繪角度切入，引入傳統(tǒng)光學(xué)攝影測(cè)量的外方位元素，建立了外方位元素整體誤差模型。發(fā)現(xiàn)3個(gè)線元素誤差導(dǎo)致的相位中心偏移在一個(gè)線性陣列的不同陣元之間是一樣的，可以借助傳統(tǒng)的SAR運(yùn)動(dòng)補(bǔ)償方法來(lái)實(shí)現(xiàn)[27]。3個(gè)角元素誤差導(dǎo)致的相位中心偏移在陣元間則不同，偏角誤差不會(huì)影響天線陣元到地面點(diǎn)的距離歷程，對(duì)成像誤差沒(méi)有影響[23]；旋角誤差會(huì)同時(shí)引入沿航向和跨航向兩個(gè)方向的運(yùn)動(dòng)誤差，補(bǔ)償較為復(fù)雜；傾角誤差并不會(huì)引入沿航向的運(yùn)動(dòng)誤差，且在跨航向-高程向平面內(nèi)的成像結(jié)果聚焦正確，僅僅位置發(fā)生了旋轉(zhuǎn)。進(jìn)一步假設(shè)除傾角外的其他外方位元素均為0時(shí)，推導(dǎo)了傾角誤差模型，提出了基于傾角誤差的MIMO下視陣列SAR誤差補(bǔ)償方法，最后通過(guò)仿真試驗(yàn)驗(yàn)證了方法的有效性。

1 外方位元素整體誤差模型

為了保障陣列天線SAR的工程應(yīng)用，引入POS所能獲得的真實(shí)數(shù)據(jù)形式來(lái)分析下視陣列SAR誤差補(bǔ)償問(wèn)題。如圖1所示，POS中的IMU(inertial measurement unit，慣性測(cè)量單元)坐標(biāo)系為原點(diǎn)在IMU中心，X軸指向航向，Y軸指向飛機(jī)右側(cè)，Z軸向下的右手直角坐標(biāo)系。輸出的POS數(shù)據(jù)是基于用戶定義的參考坐標(biāo)系，通常IMU坐標(biāo)系就是默認(rèn)的參考坐標(biāo)系。

圖1 IMU外形Fig.1 IMU appearance

IMU系統(tǒng)所獲得的3個(gè)角度(偏航角、俯仰角、橫滾角)是參考坐標(biāo)系(Ob-X′Y′Z′)相對(duì)于北東地坐標(biāo)系(On-XYZ)的旋轉(zhuǎn)角。北東地坐標(biāo)系先繞Z軸順時(shí)針旋轉(zhuǎn)偏航角κ，再繞已轉(zhuǎn)了κ角的Yκ軸順時(shí)針旋轉(zhuǎn)俯仰角ω，最后繞已轉(zhuǎn)了κ角和ω角Xκω軸順時(shí)針旋轉(zhuǎn)橫滾角α，得到的坐標(biāo)軸方向與參考坐標(biāo)系完全相同。其中，角度符號(hào)定義滿足右手規(guī)則：右手握住旋轉(zhuǎn)軸，大拇指指向旋轉(zhuǎn)軸正向，其余4指的指向就是旋轉(zhuǎn)角度的正方向。

由上述轉(zhuǎn)換關(guān)系可得北東地坐標(biāo)系與參考坐標(biāo)系之間的關(guān)系

(1)

(2)

對(duì)等式(1)兩邊求逆可得

(3)

矩陣Rα、Rω、Rκ是正交陣，因此有

(4)

通過(guò)式(4)可發(fā)現(xiàn)，該轉(zhuǎn)換關(guān)系類似于光學(xué)攝影測(cè)量中外方位元素“以X軸為主軸的轉(zhuǎn)角系統(tǒng)”。因此，為了能夠充分利用POS所獲得的數(shù)據(jù)信息，可以借鑒光學(xué)攝影測(cè)量中外方位元素的思想，為陣元位置引入3個(gè)線元素誤差和3個(gè)角元素誤差。

圖2 外方位元素誤差示意圖Fig.2 Sketch map of elements of exterior orientation errors

由IMU系統(tǒng)的旋轉(zhuǎn)角關(guān)系可知，3個(gè)角元素誤差引入的旋轉(zhuǎn)矩陣為

(5)

借鑒光學(xué)攝影測(cè)量中“空間相似變換公式”，除去比例縮放因子后可得發(fā)射陣元的實(shí)際坐標(biāo)

(6)

式中，v為載機(jī)速度；ta為方位向慢時(shí)刻；PRF為脈沖重復(fù)頻率；ΔX、ΔY和ΔZ分別為方位向、跨航向和高程向的坐標(biāo)增量。

同理可得接收陣元的實(shí)際坐標(biāo)

(7)

式中

(8)

理想情況下，等效相位中心對(duì)應(yīng)目標(biāo)點(diǎn)的距離歷程R為

(9)

(10)

根據(jù)泰勒近似公式，將兩個(gè)距離分別在慢時(shí)間ta0=x·PRF/v(x為方位向坐標(biāo))處展開可得

(11)

(12)

式中

(13)

進(jìn)一步聯(lián)立式(11)和式(12)可計(jì)算得到回波相位誤差ΔΨ為

(14)

式中，λ為波長(zhǎng)。

由式(14)可知，回波信號(hào)中的相位誤差相當(dāng)復(fù)雜，不僅與外方位元素的6個(gè)量有關(guān)，還具有3個(gè)方向的空變性，具有嚴(yán)重的耦合現(xiàn)象，需進(jìn)行分解，逐個(gè)研究單元素誤差補(bǔ)償方法。

2 傾角誤差模型

為了方便分析，假設(shè)下視陣列SAR系統(tǒng)只存在傾角誤差αy，如圖3所示?？绾较蚍植嫉膎個(gè)線陣陣元經(jīng)過(guò)等效相位中心誤差補(bǔ)償后可以等效為陣元間隔均勻的等效虛擬陣列，等效相位中心坐標(biāo)為(x,yn,0)，目標(biāo)點(diǎn)P的坐標(biāo)為(x0,y0,z0)，坐標(biāo)軸向仍然沿用圖2表述。

圖3 MIMO下視陣列SAR傾角誤差模型Fig.3 Inclination angle error model of MIMO downward-looking array SAR

借鑒光學(xué)攝影測(cè)量外方位角元素旋轉(zhuǎn)矩陣的思想，將平臺(tái)坐標(biāo)系繞X軸旋轉(zhuǎn)，則旋轉(zhuǎn)前后等效相位中心的坐標(biāo)有如下關(guān)系

(15)

因此旋轉(zhuǎn)αy角后，等效相位中心的坐標(biāo)(x′,y′,z′)變?yōu)?x,yncosαy,ynsinαy)。理想情況下，等效相位中心到目標(biāo)點(diǎn)的距離R為

(16)

(17)

則由傾角誤差引起的距離誤差ΔR為

(18)

式中

2z0ynsinαy=2y0yn(cosαy-1)+2z0ynsinαy

(19)

代入式(18)可得

(20)

由式(20)可知，距離誤差ΔR并未在x方向上引入誤差，但受目標(biāo)坐標(biāo)y0和z0的空變影響，不能直接對(duì)全孔徑進(jìn)行補(bǔ)償處理，需要考慮子孔徑補(bǔ)償方法。經(jīng)解調(diào)和匹配濾波后，回波Echo(t,x,yn)如式(21)所示。

Echo(t,x,yn)=sinc[t-2(R+ΔR)/c]exp(jΨ)=

sinc[t-2(R+ΔR)/c]

exp(-j4π(R+ΔR)/λ)

(21)

式中，sinc為sinc函數(shù)；Ψ為回波相位誤差。

3 傾角誤差補(bǔ)償方法

MIMO下視陣列SAR地形重建最直接的方法是利用等效相位中心原理，將多發(fā)多收收發(fā)分置的回波信號(hào)轉(zhuǎn)化為收發(fā)共用的情況后，再沿用傳統(tǒng)SAR二維成像技術(shù)和波束形成技術(shù)實(shí)現(xiàn)三維成像。當(dāng)載機(jī)平臺(tái)受外界環(huán)境等影響產(chǎn)生運(yùn)動(dòng)誤差后，成像補(bǔ)償方法習(xí)慣將該誤差利用三角函數(shù)轉(zhuǎn)換為沿航跡方向和雷達(dá)視線方向位置位移來(lái)分別處理。沿航跡方向非均勻采樣問(wèn)題通常用時(shí)域插值重采樣方法解決；雷達(dá)視線方向補(bǔ)償則多忽略目標(biāo)空變性影響，直接進(jìn)行全孔徑補(bǔ)償處理。首先采用場(chǎng)景中心線作為參考目標(biāo)對(duì)整個(gè)場(chǎng)景進(jìn)行補(bǔ)償，然后對(duì)參考線以外的其他距離目標(biāo)回波的殘留相位誤差進(jìn)行補(bǔ)償。但從傾角誤差模型來(lái)看，距離誤差ΔR受目標(biāo)坐標(biāo)y0和z0的空變影響，不能直接進(jìn)行全孔徑補(bǔ)償。

為了消除距離誤差的空變影響，將式(21)中的相位分別對(duì)x和yn求偏導(dǎo)，得到航跡向波數(shù)kx和跨航向波數(shù)ky

(22)

(23)

聯(lián)立式(22)、式(23)和式(17)

(24)

利用菲涅爾定理化簡(jiǎn)，并忽略陣元長(zhǎng)度小值，可得

(25)

代入式(24)可得

(26)

將式(25)和式(26)代入式(20)可得距離誤差

(27)

由式(27)可知，距離誤差與目標(biāo)位置無(wú)關(guān)，即消除了目標(biāo)空變影響。將其分為ΔR1和ΔR2兩部分，ΔR1與波數(shù)無(wú)關(guān)，可在回波的三維空域中首先對(duì)其進(jìn)行補(bǔ)償，稱為一次補(bǔ)償。ΔR2與波數(shù)相關(guān)，可在二維波數(shù)域中分塊計(jì)算，稱為二次補(bǔ)償。

首先由ΔR1構(gòu)造濾波器Hα1，進(jìn)行一次補(bǔ)償。

Hα1=exp(-j4πΔR1/λ)=

(28)

再將補(bǔ)償完ΔR1的回波數(shù)據(jù)做航跡向和跨航向傅里葉變換，變換到該二維波數(shù)域，并沿航跡向和跨航向均勻分為M和N段，獲得M×N個(gè)子塊，每個(gè)子塊用中心波數(shù)kx-ic和ky-jc代替kx和ky計(jì)算距離誤差ΔR2-ij為

(29)

將每個(gè)子塊補(bǔ)零并擴(kuò)展到原尺寸后，再做航跡向和跨航向傅里葉逆變換，以及距離向傅里葉變換，得到航跡向和跨航向二維空域，距離向波數(shù)域回波

(30)

式中，kr為距離向波數(shù)；krc為中心波數(shù)；Br為帶寬。

將式(30)乘以濾波器Hα2來(lái)補(bǔ)償ΔR2-ij引起的相位誤差，稱為二次補(bǔ)償。再做距離向傅里葉逆變換，轉(zhuǎn)換到三維空域，把每個(gè)子塊補(bǔ)償?shù)臄?shù)據(jù)直接相加，得到最終回波。接著使用三維RD算法進(jìn)行成像即可得到理想的成像結(jié)果。傾角誤差補(bǔ)償及成像流程如圖4所示。

(31)

圖4 傾角誤差補(bǔ)償及成像流程Fig.4 Flow of inclination angle error compensation and imaging

4 試驗(yàn)與分析

4.1 仿真地形試驗(yàn)

本文仿真了一個(gè)具有3個(gè)山峰、兩個(gè)山谷的山區(qū)地形，采用2發(fā)80收的MIMO模式，利用表1中的仿真參數(shù)對(duì)具有傾角誤差的MIMO下視陣列SAR系統(tǒng)進(jìn)行三維成像試驗(yàn)。圖5(a)和5(b)分別給出原始仿真地形的三維圖和二維圖，最大高程為48.626 0 m，最小高程為-39.297 7 m。

表1 載機(jī)系統(tǒng)仿真參數(shù)

圖5 原始仿真地形Fig.5 Original simulation terrain

航測(cè)規(guī)范中要求傾角一般不大于2°，最大不超過(guò)3°，因此仿真試驗(yàn)時(shí)若故意把傾角加大來(lái)分析，最終結(jié)果并沒(méi)有意義。本文在仿真過(guò)程中選取最大傾角為3°，如圖6所示。理想情況下，在YOZ平面陣元的位置為(nd,0)(n代表第n個(gè)等效虛擬陣元，d代表等效虛擬陣元之間的間隔)，當(dāng)載機(jī)飛行過(guò)程中發(fā)生了傾角誤差時(shí)，陣元的位置在Y方向和Z方向則會(huì)分別產(chǎn)生分量ndcosαy和ndsinαy。

圖6 傾角誤差帶來(lái)陣元位置變化示意Fig.6 Sketch map of array element position change caused by inclination angle error

補(bǔ)償前，陣列SAR原始回波數(shù)據(jù)如圖7(a)所示。成像處理后實(shí)現(xiàn)斜距、方位和角度三維分辨，圖7(b)、圖7(c)、圖7(d)分別給出方位-角度平面、角度-斜距平面和方位-斜距平面的切面圖，從目視效果來(lái)看，與角度相關(guān)的切面圖均出現(xiàn)明顯的整體位移。將成像結(jié)果根據(jù)角度、斜距與高度、地距之間的轉(zhuǎn)換關(guān)系，按照與原始地形相同的格網(wǎng)間隔插值規(guī)劃到同一坐標(biāo)系后，可得到三維成像結(jié)果圖7(e)，很明顯仿真區(qū)域整體產(chǎn)生了偏移，且進(jìn)一步導(dǎo)致觀測(cè)區(qū)域邊緣產(chǎn)生較大誤差。圖7(f)給出了成像結(jié)果的二維圖，與圖5(b)相比，整體效果欠佳，由-39.30至48.63的高程范圍增大到-39.00至77.13之間，誤差較大。

經(jīng)等效相位中心相位誤差補(bǔ)償后的回波數(shù)據(jù)如圖8(a)所示，振幅大小有了明顯地縮小，有效補(bǔ)償了相應(yīng)回波。經(jīng)過(guò)成像處理后，圖8(b)、圖8(c)、圖8(d)分別給出補(bǔ)償后方位-角度平面、角度-斜距平面、方位-斜距平面的切面圖。方位-角度平面由圖7(b)中跨航向的15 m～125 m區(qū)域校正到25 m～135 m區(qū)域，角度-斜距平面由7(c)中沿航向的10 m～120 m區(qū)域校正到30 m～140 m區(qū)域。將成像結(jié)果按相同格網(wǎng)間隔插值規(guī)劃到地距、方位、高度坐標(biāo)系后，得到圖8(e)所示的三維成像結(jié)果，補(bǔ)償前整幅圖像上出現(xiàn)的大量毛刺信息消失了。與圖7(f)相比，成像結(jié)果二維圖8(f)更能明顯地展示出改善效果，3個(gè)山峰和2個(gè)山谷都基本恢復(fù)了本來(lái)形貌。

圖7 仿真地形補(bǔ)償前結(jié)果圖Fig.7 Pre-compensation result maps of simulation terrain

圖8 仿真地形補(bǔ)償后結(jié)果圖Fig.8 Compensated result maps of simulation terrain

為了進(jìn)一步定量分析補(bǔ)償效果，將補(bǔ)償前的成像結(jié)果(圖7(e))與原始地形(圖5(a))作差分處理后可得差分圖9(a)，在山峰和山谷地區(qū)出現(xiàn)多個(gè)明顯的細(xì)胞狀誤差區(qū)域；統(tǒng)計(jì)其直方圖10(a)，差值在-31.051 m至25.121 m之間變換，且有大量不在0附近的值，誤差較大。將補(bǔ)償后的成像結(jié)果(圖8(e))與原始地形(圖5(a))作差分處理后可得差分圖9(b)，很明顯高程誤差得到了較好地控制；統(tǒng)計(jì)其直方圖10(b)，差值縮小至[-4.356 m，4.205 m]區(qū)間，且主要集中在-1 m至1 m之間。但從整體上來(lái)看，圖面信息仍較亂，尤其在幾個(gè)山峰和山谷地區(qū)，有大量繞著山勢(shì)類似等高線的誤差區(qū)域。究其原因，除了跨航向分辨率的局限以外，主要還是距離計(jì)算不準(zhǔn)確所致，后續(xù)可考慮借助控制信息校正相應(yīng)距離值。

圖9 仿真地形的成像結(jié)果與原始場(chǎng)景差分圖Fig.9 Difference map between imaging results and original scene of simulation terrain

圖10 差分結(jié)果直方圖Fig.10 Differential results histogram

表2定量統(tǒng)計(jì)了補(bǔ)償前后成像結(jié)果的高程誤差均值、標(biāo)準(zhǔn)差和整個(gè)場(chǎng)景高程誤差在半個(gè)分辨率之內(nèi)的概率。補(bǔ)償后，前兩個(gè)指標(biāo)有了明顯下降，最后一個(gè)指標(biāo)又有大幅提高，這足以說(shuō)明補(bǔ)償方法的有效性。

表2 成像結(jié)果定量指標(biāo)統(tǒng)計(jì)

4.2 SRTM DEM真實(shí)地形試驗(yàn)

為了進(jìn)一步驗(yàn)證MIMO下視陣列SAR傾角誤差補(bǔ)償方法的有效性，本文選取青藏高原地區(qū)30 m分辨率的SRTM DEM為原始場(chǎng)景，仍然采用2發(fā)80收的MIMO模式進(jìn)行試驗(yàn)。圖11給出對(duì)原始DEM重采樣為0.5 m格網(wǎng)間距后截取的200×200 m大小的DEM二維圖和三維圖。最大高程為4 931.2 m，最小高程為4 712.1 m。由于地形高程的抬高，表1中平臺(tái)高度也相應(yīng)調(diào)整為5200 m，其他參數(shù)不變。

圖11 SRTM DEM真實(shí)地形Fig.11 SRTM DEM real terrain

補(bǔ)償前的三維成像結(jié)果如圖12所示，不僅整體位置發(fā)生偏移；大量高程值也不準(zhǔn)確，尤其在跨航向上發(fā)生偏移的反方向表現(xiàn)明顯。利用本文方法補(bǔ)償后，三維成像結(jié)果如圖13所示，基本恢復(fù)了真實(shí)DEM的本來(lái)形貌。為了能更加直觀地表現(xiàn)補(bǔ)償效果，分別將補(bǔ)償前后的成像結(jié)果與原始SRTM DEM作差分，可得圖14，誤差范圍從[-151.56，133.15]縮小到[-9.99，9.09]。統(tǒng)計(jì)表2中的3項(xiàng)指標(biāo)，高程誤差均值從-24.05降到-0.03；高程誤差標(biāo)準(zhǔn)差從49.88降到1.82；整個(gè)場(chǎng)景高程誤差在半個(gè)分辨率之內(nèi)的概率從2.38%提高到80.28%。無(wú)論從目視效果，還是統(tǒng)計(jì)量都充分證明了補(bǔ)償方法的有效性。

5 結(jié)束語(yǔ)

陣列天線子陣相位中心高精度測(cè)量是MIMO下視陣列SAR系統(tǒng)的重點(diǎn)。系統(tǒng)平臺(tái)實(shí)際運(yùn)動(dòng)中產(chǎn)生的位置偏移和姿態(tài)變化會(huì)造成天線相位中心的偏移，影響高分辨率陣列SAR三維圖像質(zhì)量。本文結(jié)合POS真實(shí)數(shù)據(jù)形式，引入傳統(tǒng)光學(xué)攝影測(cè)量的外方位元素，提出了基于傾角誤差的MIMO下視陣列SAR誤差補(bǔ)償方法。該方法借助回波相位在沿航向和跨航向的偏導(dǎo)值，以及存在傾角誤差時(shí)的斜距值聯(lián)立求解，消除了傾角誤差的空變影響，采用二維波數(shù)域分塊計(jì)算距離誤差，空域逐塊補(bǔ)償?shù)姆绞酵瓿闪巳S成像。通過(guò)仿真試驗(yàn)驗(yàn)證了本文方法的有效性。

后續(xù)研究中，需進(jìn)一步評(píng)定誤差補(bǔ)償后的精度，通過(guò)對(duì)陣列SAR數(shù)據(jù)處理中若干非理想因素進(jìn)行補(bǔ)償和影響分析，明確指出對(duì)載機(jī)POS的要求。

圖12 真實(shí)地形補(bǔ)償前結(jié)果圖Fig.12 Pre-compensation result maps of real terrain

圖13 真實(shí)地形補(bǔ)償后結(jié)果圖Fig.13 Compensated result maps of real terrain

圖14 真實(shí)地形的成像結(jié)果與原始DEM差分圖Fig.14 Difference map between imaging results and original DEM of real terrain