師 君 馬 龍 韋順軍 時代奇 張曉玲 陳 剛

①(電子科技大學電子工程學院 成都 611731)

②(西安測繪研究所 西安 710054)

基于導航數(shù)據(jù)的Ka波段InSAR成像處理與分析

師 君*①馬 龍①韋順軍①時代奇①張曉玲①陳 剛②

①(電子科技大學電子工程學院 成都 611731)

②(西安測繪研究所 西安 710054)

與其它波段相比，毫米波干涉SAR (InSAR)具有精度高、體積小的特點，是目前干涉SAR研究的重要方向之一。但由于較短的波長，其也給成像處理和干涉相位提取帶來了一定困難。該文利用毫米波段干涉SAR數(shù)據(jù)，采用流結構后向投影(BP)算法結合IMU進行了毫米波SAR成像和干涉相位分析。研究發(fā)現(xiàn)，由于毫米波波長短，對天線相位中心變化更為敏感，為了保證系統(tǒng)的干涉相位，必須利用IMU數(shù)據(jù)進行精確的運動誤差補償。通過實測數(shù)據(jù)驗證了采用各自計算天線相位中心的策略，后向投影算法可以在成像過程中精確補償平地相位。

干涉SAR成像；成像性能分析；后向投影算法；基于導航數(shù)據(jù)的運動補償

1 引言

與其它波段相比，Ka波段SAR穿透能力弱，能更好地反應積雪、植被等地貌的高程信息，且系統(tǒng)結構緊湊，重量輕。因此，近年來國外紛紛開展了Ka波段干涉SAR系統(tǒng)研制[1－5]。但另一方面，由于Ka波段波長更短，為了獲得良好的干涉相位，需要系統(tǒng)具有較好的通道一致性且成像方法具有較好的相位保持能力。

對于機載 SAR系統(tǒng)，由于風場、湍流等因素的影響，天線相位中心軌跡較理想軌跡可能存在較大的偏差，嚴重影響SAR(特別是高波段SAR)圖像質(zhì)量，進而影響高精度DEM提取等SAR應用。由于傳統(tǒng)的距離-多普勒(RD)算法、CS算法都基于SAR的多普勒質(zhì)心、多普勒斜率模型，當IMU測量數(shù)據(jù)較為復雜時，上述方法需要較大改進才能實現(xiàn)對復雜運動誤差的精確補償。與之相比，后向投影(BP)算法[6－11]采用逐點匹配的處理方法，可實現(xiàn)IMU數(shù)據(jù)與成像處理的充分結合，更適合于復雜運動SAR成像處理。

本文基于實測Ka波段毫米波干涉SAR數(shù)據(jù)，對后向投影算法的在高分辨率 SAR成像以及干涉相位提取方面的性能進行了分析。第2節(jié)首先回顧了流結構BP算法的處理流程，并針對所處理數(shù)據(jù)由于距離向去斜率處理產(chǎn)生的存在的殘留相位問題，對BP算法的補償相位進行了修正。第3節(jié)從成像的角度出發(fā)，分析了實測數(shù)據(jù)的回波特征，IMU數(shù)據(jù)對成像結果的影響，以及去斜率處理對結果的影響。第4節(jié)從干涉相位的角度出發(fā)，對IMU數(shù)據(jù)、去斜率處理等對BP算法的保相性能進行了分析，并初步驗證了地距投影條件下，BP算法所特有的去平地能力。

2 流結構后向投影算法

2.1 算法介紹

后向投影算法(BP)通過將回波數(shù)據(jù)逐點投影到圖像空間各像素，實現(xiàn)各散射點能量的積累。與其它方法相比，BP方法成像原理簡單，處理過程不存在任何近似誤差，可適用于各種模式 SAR成像處理，且便于進行高精度運動誤差補償。BP算法工作示意圖如圖1所示。

一般地，BP算法包括如下步驟：

步驟1 確定圖像空間

可以選擇距離-方位向作為投影空間，以使得BP算法成像結果與其它算法成像結果兼容，也可以選擇地面作為投影空間，以便于后續(xù)圖像處理，像素間隔應略小于系統(tǒng)理論分辨率。

步驟2 距離壓縮與插值

為了提高成像質(zhì)量，需要對距離壓縮后數(shù)據(jù)進行插值，可選擇頻域sinc插值。

步驟3 計算距離歷史

給定某像素p，其距離歷史為：

其中，pAPC(n)表示天線相位中心在第n個PRI的位置，表示向量2范數(shù)，n為慢時間域，p為散射點位置。

步驟4 計算回波位置

根據(jù)圖1所示的雷達回波時序，某像素在第n個PRI對應的回波位置為：

其中，[x]表示4舍5入取整，TDelay表示接收波門延時，c為光速，TPulse表示發(fā)射脈沖持續(xù)時間，κ表示頻域插值倍數(shù)，fs表示采樣頻率。

步驟5 相參積累

補償該散射點的多普勒相位，并將其投影到圖像空間中，

其中，fc為發(fā)射信號載頻，I(x,y)表示圖像空間，τ(n,p)為散射點到雷達的延時，“←”表示賦值操作。

步驟6 遍歷所有像素和PRI

通過上述操作，即可實現(xiàn)對處理區(qū)域的成像處理。

2.2 去斜處理改進

由于目前實測數(shù)據(jù)距離壓縮采用了去斜率處理技術，而去斜率處理會導致額外的方位向相位，因此在BP成像過程中，需要對殘留相位進行補償。假設系統(tǒng)發(fā)射線性調(diào)頻信號，則單點目標回波可表示為：

其中，t為快時間域，K為發(fā)射信號的調(diào)頻斜率，

去斜處理將回波信號與發(fā)射信號共軛相乘，即：

將線性調(diào)頻信號相位的平方項展開，可以得到去載頻和去斜率后的信號為：

圖1 BP算法工作示意圖Fig. 1 Illustration on BP algorithm

式(7)中包含 3個復指數(shù)分量。第 1項為傳統(tǒng)SAR系統(tǒng)的多普勒相位，產(chǎn)生方位向分辨率。第2項為快時間t的單頻信號，且頻率與目標的延時相關，因此，通過傅里葉變換即可實現(xiàn)對去斜率處理信號的距離壓縮以及不同距離目標的分辨。去斜率處理距離向的分辨率由傅里葉變換的頻率分辨率決定，并與發(fā)射信號的調(diào)頻斜率和采樣時間有關，兩者乘積即為發(fā)射信號的帶寬。因此，去斜率處理的距離向分辨率與脈沖壓縮系統(tǒng)相同，ρR=c/(2B)。

式(7)中的第 3項為去斜率處理產(chǎn)生的殘留相位，由于該項與快時間t無關，因此該項不會對距離壓縮產(chǎn)生影響。但是另一方面，由于其與慢時間n有關，因此該項將會對系統(tǒng)的方位向處理產(chǎn)生影響，需要在相位補償過程中加以補償。此時，后向投影算法的相位補償公式為：

3 成像性能分析

為了獲得較高的干涉相位，需要首先獲得高精度的各通道 SAR圖像，本節(jié)主要從信號特征，運動補償和殘留相位3個方面，分析后向投影算法的成像性能。

3.1信號分析

為了分析圖像中強點目標的信號特征，首先根據(jù)天線相位中心軌跡，采用式(2)計算強點在距離壓縮后數(shù)據(jù)中的位置，并逐一提取各位置處的復數(shù)據(jù)，如圖2(a)所示，可以看出，計算所得ID(紅色線段)與實測的強點回波良好吻合。

圖2(b)為強點目標的方位向頻譜，可以看出其表現(xiàn)出了 SAR方位向信號典型的帶限信號特征。圖2(c)為強點目標的時頻譜，可以看出其表現(xiàn)出了SAR方位向信號典型的線性調(diào)頻信號特征。圖2(d)為強點目標的DCFT譜(Dechirp變換)，可以看出其與理想線性調(diào)頻信號的DCFT譜類似，但并非精確的十字結構，表明該強點可能并非嚴格意義上的單點目標。

3.2 IMU數(shù)據(jù)分析

假設 SAR系統(tǒng)天線相位中心實際軌跡為，給定散射點p，回波對應的實際距離歷史為：

圖2 Ka波段SAR強點目標的回波特征Fig. 2 Signal characteristics of Ka-band SAR data

在實際軌跡未知的條件下，成像過程往往假設平臺直線運動。此時，用于補償?shù)木嚯x歷史為：

和之間的誤差即為運動誤差，其導致兩方面的影響：在補償距離走動效應時，計算ID可能存在偏差，導致遺漏散射點的部分回波，并進入了額外的噪聲數(shù)據(jù)；在計算補償相位時出現(xiàn)誤差，導致無法完全補償散射點的多普勒相位，影響圖像聚焦和干涉相位精度。其中，第1類誤差對距離歷史精度要求為距離分辨率級，一般為分米級；第2類誤差對距離歷史精度要求為亞波長級，一般需要達到毫米級。

為了進一步分析，可以對式(9)和式(10)做差，得到距離歷史誤差近似公式為：

其中，e(n)為天線相位中心實際軌跡與直線近似軌跡之差。

式(11)的第1項為平臺到散射點的方向矢量，式(11)的物理意義十分清晰，即距離誤差為天線相位中心軌跡誤差沿距離向的投影。按照合成孔徑雷達理論，為了實現(xiàn)精確成像，相位誤差應控制在波長的1/8以內(nèi)，對于Ka波段，距離誤差應小于0.5 mm。圖3為IMU測量得到的實際誤差?？梢钥闯?，平臺較直線軌跡的偏移量達到了至少分米量級。

因此，必須利用 IMU測量數(shù)據(jù)進行補償。相應地，根據(jù)SAR相位補償對誤差要求，IMU測量精度也應達到毫米量級。需要說明的是，1/8的相位誤差條件基于隨機噪聲相位模型。針對 IMU測量系統(tǒng)，需要保證預處理后(濾波、擬合等)的天線相位中心軌跡足夠光滑，隨機抖動小于1 mm。對于 IMU測量單元的隨機游走類偏差以及卡爾曼濾波導致的偏差，仍需要進一步根據(jù)誤差類型加以分析。保守地說，上述誤差在一個孔徑時間內(nèi)的積累應小于1 mm(此時，上述誤差可以完全忽略不計，但誤差如果大于1 mm，還需要根據(jù)誤差的時間特性(線性特征、高階特征、隨機特征)加以進一步分析與討論。

圖4(a)和圖5(a)為采用平均速度進行BP成像的結果(載頻約為35 GHz，信號帶寬約1 GHz，平臺速度約60 m/s，作用距離約1200 m，入射角約45°，場景大小8192像素×8192像素，距離向像素間隔0.10 m，方位向像素間隔0.05 m)，可以看出，雖然圖像實現(xiàn)的聚焦，但是整幅圖像的對比度較低，強點目標在圖像中存在較為明顯的散焦，且圖中的道路等目標輪廓較為模糊。

圖3 航跡高階誤差Fig. 3 High-order trajectory errors

圖4 整體效果對比Fig. 4 Comparison on the whole image

圖4(b)和圖5(b)為利用IMU數(shù)據(jù)進行非勻速運動誤差補償后的成像結果，可以看出，此時圖像的對比度更大，局部區(qū)域的道路輪廓更為明顯，且強點目標散焦現(xiàn)象較圖5(b)大大減弱，但存在一定程度的散焦。

3.3 去斜率影響分析

由于目前階段的系統(tǒng)在距離向采用了去斜率處理技術，導致圖像在方位向存在一定的殘留相位，可能會影響最終成像質(zhì)量。因此，本節(jié)對殘留相位的影響進行了分析。

圖6為利用式(7)得到的多普勒分量和殘留相位分量產(chǎn)生的方位向信號頻譜。從中可以看出，多普勒分量導致的方位向頻譜有大約1200 Hz，而由于殘留相位產(chǎn)生的方位向頻譜只有大約3 Hz左右。因此，多普勒分量在目前系統(tǒng)的方位向分辨率中占主要貢獻，殘留相位的影響較小。

圖7為補償殘留相位后得到的BP成像結果，可以看出，補償殘留相位對圖像整體效果影響不大。但對比局部圖可以發(fā)現(xiàn)，補償殘留相位后圖像的對比度要略微優(yōu)于未補償圖像。

圖8為補償殘留相位前后某強點目標成像效果比較，可以看出，目前成像處理對于強點目標仍存在一定的殘留距離走動和方位向散焦，其可能原因首先是目前圖像像素間隔很密(0.10 m×0.05 m)，對系統(tǒng)參數(shù)、運動誤差補償?shù)囊筝^高，可能存在參數(shù)估計不準的問題。對比補償殘留相位前后的圖像，可以發(fā)現(xiàn)，未補償殘留相位時，該強點散焦能量向左側偏移，而補償殘留相位后，強點散焦能量較為對稱。因此，在干涉處理中目前選擇補償殘留相位后圖像進行分析(研究過程中對兩種補償方法的干涉相位進行了分析，發(fā)現(xiàn)是否補償該相位對兩幅圖像的相參性及殘差點數(shù)目影響不大)。

4 相位性能分析

本節(jié)利用Ka波段干涉SAR的兩幅圖像對目前BP成像算法的相位特性進行分析。

4.1 BP算法的去平地能力

圖5 局部區(qū)域?qū)Ρ菷ig. 5 Comparison on the local region

圖6多普勒相位與殘留相位方位向頻譜Fig. 6 Frequency spectrums of the Doppler term and residualterm

傳統(tǒng)的SAR成像算法基于多普勒質(zhì)心和多普勒斜率模型，即通過補償距離歷史的1階分量、2階分量及高階分量實現(xiàn)對散射點的能量聚集。由于實際距離歷史的 0階分量被保留，因此干涉相位為主、輔天線到散射點最短距離(理想正側視條件下)之差。

圖7 補償殘留相位后BP成像結果Fig. 7 Imaging result after compensating the residual phase

圖8 補償殘留相位前后成像結果對比Fig. 8 Comparison on the imaging results with and without the residual phase

而對于BP算法，由于采用精確距離歷史進行相位補償，因此，用于進行干涉處理的0階分量也在成像過程中剔除。為了保證BP算法與其它算法兼容，最簡單的方法是在補償干涉相位時保留距離歷史的0階相位、此時，補償相位為(此處忽略由于去斜率引入的殘留相位，并將相位表示為波數(shù)的形式)：

其中，K0為載波波數(shù)，R(p,0)表示中心時刻的斜距。

當忽略平臺的姿態(tài)誤差時，主天線的距離歷史可以近似認為與輔天線的距離歷史平行。此時，主天線BP與輔天線BP具有相同的多普勒斜率及高階項。采用標準的BP算法(不減去0階項)，利用主天線的參數(shù)可以實現(xiàn)輔天線通道成像處理。

觀察式(12)可以發(fā)現(xiàn)，由于在成像過程中都減去了主天線的0階分量(注意，如果主輔天線利用各自參數(shù)進行成像處理，則是減去各自天線的0階分量)。該過程相當于利用式(12)進行成像處理后，在主輔天線上各乘了一個相同的相位。由于干涉相位提取時需要對主輔圖像相位做差，相同的常數(shù)相位不會對成像過程產(chǎn)生影響。

因此，忽略主輔天線相位中心的位置差別，只利用主天線的參數(shù)對兩幅圖像進行成像處理得到的干涉相位與其它成像處理算法得到的干涉相位相同。同樣地，此時BP算法也存在由于基線因素導致的圖像偏移，需要在干涉相位提取中采用圖像配準、去平地等技術進行補償。

為了進一步提高成像處理精度，需要在處理過程中考慮姿態(tài)誤差，此時主輔天線相位中心軌跡不再平行，無法統(tǒng)一利用主天線進行成像處理，而需要利用主輔天線各自參數(shù)分別成像。

為了得到與傳統(tǒng)干涉 SAR一致的相位，則需要補償?shù)闹鬏o天線相位項分別為：

其中，M,S分別代表主天線和輔天線。如果不補償0階項，則干涉相位為：

其中，?In表示式(13)對應的理想的干涉SAR相位，?BP表示標準 BP算法分別用主輔天線參數(shù)成像得到的干涉相位。

觀察最后一項可以發(fā)現(xiàn)，?In和?BP之間相差的相位為假設散射點位于像素點時(一般情況下，散射點位于3維地理空間，像素點位于2維圖像空間，散射點被“映射”到像素點，而非實際“位于”像素點)的干涉相位。由于像素點的高度為 0(選擇水平面作為BP算法投影面)，所以該相位實際上是像素點位置的平地相位。

因此，通過在成像處理過程中精確計算主天線和輔天線的天線相位中心歷史，并利用各自參數(shù)進行成像處理，可以剔除干涉相位中的平地相位。另外，BP算法的投影特性還可以補償部分的主輔圖像間的位置偏差，但無法補償由于高程導致的主副圖像像素偏移(由于散射點高度未知)。

4.2 實測數(shù)據(jù)分析

圖 9為采用上述兩種策略得到的干涉相位結果比較。圖9(a)為只利用主天線進行雙通道數(shù)據(jù)成像得到的干涉相位圖(兩幅圖的相參性為 0.9174 (未剔除遮擋區(qū)域像素點)，殘差占整幅圖比例8.865%)，可以看出，此時圖像存在明顯的平地相位，需要進行相應的去平地處理。圖9(b)為精確計算各自天線相位中心，并進行補償，得到的干涉相位圖(兩幅圖的相參性為 0.9156，殘差占整幅圖比例9.110%)?？梢钥闯?，此時干涉圖中的平地相位并不顯著，可以直接進行濾波、解纏等后續(xù)處理。圖9(c)為圖9(b)濾波后的干涉相位圖，能反映出典型的地形特征。

圖 10為圖 9(c)局部區(qū)域的干涉相位圖和對應的3維效果圖(目前處理未包含相位解纏操作)，從中可以看出，對于A區(qū)域，干涉相位圖中存在顯著的條狀紋理特征，能反映出實際地面上的規(guī)則高低起伏，對于B區(qū)域，干涉相位圖中存在顯著的塊狀紋理特征，也反映出實際地面上的規(guī)則高低起伏。

5 結論

通過對實測毫米波InSAR數(shù)據(jù)分析，可以看出：

(1) 目前的毫米波InSAR具有較好的相參性，能夠得到較高精度的干涉相位。

(2) 由于毫米波波長短，對天線相位中心變化更為敏感，為了保證系統(tǒng)的干涉相位，必須利用IMU數(shù)據(jù)進行精確的運動誤差補償。

圖9 BP算法得到的干涉相位圖(方位向隔一抽取，保證整幅圖像0.1 m×0.1 m比例)Fig. 9 Interferograms by using BP algorithm (down-sampling (base-2) in the azimuthdirection to ensure 0.1 m×0.1 m scale)

圖10 局部干涉相位圖Fig. 10 Interferograms of local regions

(3) 通過實測數(shù)據(jù)，驗證了采用各自計算天線相位中心的策略可以在成像過程中精確補償平地相位。

(4) 從干涉相位也可以看出，目前的干涉相位圖仍存在較大的噪聲，其可能由于未進行精確圖像配準、未進行精確運動誤差補償?shù)纫蛩匾?。因此需要進一步開展深入研究。

致謝感謝中航科工二院23所提供研究所需的Ka波段干涉SAR數(shù)據(jù)。

[1] D’Addio S and Ludwig M. Modeling analysis of rain effect on Ka-band single pass InSAR performance[C]. IEEE International Geoscience and Remote Sensing Symposium (IGARSS2009), Cape Town, Republic of South Africa 2009, 4: 913-916.

[2] Ludwig M, D’Addio S, and Engel K. A spaceborne Ka-band SAR interferometer concept based on scan-on-receive techniques[C]. 2010 8th European Conference on Synthetic Aperture Radar (EUSAR), Berlin, Germany, 2010: 1-4.

[3] Schmitt M and Stilla U. Adaptive multilooking of airborne Ka-band multi-baseline InSAR data of urban areas[C]. 2012 IEEE International Geoscience and Remote Sensing Symposium (IGARSS), Munich, Germany, 2012: 7401-7404.

[4] Schaefer C and Lopez-Dekker P. Interferometric Ka-band SAR with DBF capability[C]. 2012 9th European Conference on Synthetic Aperture Radar, Nuremberg, Germany, 2012: 7-10 .

[5] Mokadem A, Thirion-Lefevre L, and Colin-Koeniguer E. Analysing urban areas in the frame of non-line of sight target detection electromagnetic modelling, validation and application to real data in Ka-band[C]. 2013 IEEE International Conference on Electromagnetics in Advanced Applications (ICEAA), Torino, Italy, 2013: 543-546.

[6] Peters T M. Algorithms for fast back- and re-projection in computed tomography[J].IEEE Transactions on Nuclear Science, 1981, 28(4): 3641-3647.

[7] Ulander L M H, Hellsten H, and Stenstrom G. Synthetic-aperture radar processing using fast factorized back-projection[J].IEEE Transactions on Aerospace and Electronic Systems, 2003, 39(3): 760-776.

[8] Cantalloube H and Dubois-Fernandez P. Airborne X-band SAR imaging with 10 cm resolution: technical challenge and preliminary results[J].IEE Proceeding-Radar,Sonar and Navigation, 2006, 153(2): 163-176.

[9] Ash J N. An Autofocus method for backprojection imagery in synthetic aperture radar[J].IEEE Geoscience and Remote Sensing Letters, 2012, 9(4): 104-108.

[10] Tan Wei-xian, Li Dao-jing, and Hong Wen. Airborne spotlight SAR imaging with super high resolution based on backprojection and autofocus algorithm[C]. IEEE International Geoscience and Remote Sensing Symposium (IGARSS2008), Boston, Massachusetts, USA, 2008, 4: 1300-1303.

[11] Shi Jun, Ma Long, and Zhang Xiao-ling. Streaming BP for non-linear motion compensation SAR imaging based on GPU[J].IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2013, 6(4): 2035-2050.

師 君(1979-)，男，河南人，獲電子科技大學工學博士學位，目前為電子科技大學副教授，主要從事SAR數(shù)據(jù)處理方面研究。

E-mail: shijun@uestc.edu.cn

馬 龍(1990-)，男，安徽人，獲電子科技大學學士學位，目前在電子科技大學攻讀碩士學位，主要從事SAR成像技術研究。

E-mail: 1057616261@qq.com

韋順軍(1983-)，男，廣西人，獲電子科技大學工學博士學位，目前為電子科技大學講師，主要從事SAR成像技術、干涉SAR技術研究。

E-mail: 249785144@qq.com

時代奇(1988-)，男，河北人，獲電子科技大學學士學位，目前在電子科技大學攻讀碩士學位，主要從事干涉SAR成像技術研究。

E-mail: 654729253@qq.com

張曉玲(1964-)，女，四川人，獲電子科技大學工學博士學位，目前為電子科技大學教授/博導，主要從事 SAR成像技術、雷達探測技術研究。

E-mail: xlzhang@uestc.edu.cn

陳 剛(1976-)，男，陜西人，獲國防科技大學工學碩士學位，目前為西安測繪研究所副研究院/副主任，主要從事InSAR數(shù)據(jù)處理、干涉定標技術研究。

E-mail: splitter@263.net

Ka-band InSAR Imaging and Analysis Based on IMU Data

Shi Jun①Ma Long①Wei Shun-jun①Shi Dai-qi①Zhang Xiao-ling①Chen Gang②

①(School of Electronic Engineering, University of Electronic Science and Technology of China, Chengdu 611731, China)
②(Xi’an Research Institute of Surveying and Mapping, Xi’an 710054, China)

Compared with other bands, the millimeter wave Interferometric Synthetic Aperture Radar (InSAR) has high accuracy and small size, which is a hot topic in InSAR research. On the other hand, shorter wavelength causes difficulties in 2D imaging and interferometric phase extraction. In this study, the imaging and phase performance of the streaming Back Projection (BP) method combined with IMU data are analyzed and discussed on the basis of actual Ka-band InSAR data. It is found that because the wavelength of the Ka-band is short, it is more sensitive to the antenna phase-center history. To ensure the phase-preserving capacity, the IMU data must be used with accurate motion error compensation. Furthermore, during data processing, we verify the flat-earth-removing capacity of the BP algorithm that calculates and compensates the master and slave antenna phase centers individually.

InSAR imaging; Imaging performance analysis; Back Projection (BP) algorithm; Motion compensation based on IMU data

中國分類號：TN957.52

A

2095-283X(2014)01-0019-09

10.3724/SP.J.1300.2014.13142

2013-12-24收到，2014-02-21改回；2014-03-06網(wǎng)絡優(yōu)先出版

61101170)，博士后基金(2011018511001)和高分專項(GFZX04031102)資助課題

*通信作者： 師君 shijun@uestc.edu.cn