一種機載重軌干涉SAR相位偏移量估計算法

， ，

(1.中國科學院電子學研究所， 北京 100190；2.中國科學院大學， 北京 100049；3.空間信息處理與應用系統技術重點實驗室， 北京 100190)

0 引言

機載重軌干涉SAR由于其飛行機動靈活、圖像分辨率高、重訪周期短等優點，能彌補星載SAR固有的缺點，特別適合對局部地區進行高精度地形測繪[1-2]、形變檢測[3]和三維成像[4]，具有較高的航空遙感運用價值。

在干涉SAR反演高程的過程中，絕對相位與解纏繞后的相位存在一個常數偏移量。目前方法利用在成像場景中人工部署并測量一定數量的定標點來計算常數偏移量。然而，實際工程中人為布設定標點浪費人力財力，尤其在一些危險區域部署和測量定標點是極其困難和難以實現的。為了解決這類問題，目前自動計算相位偏置的方法一般分為基于數據冗余的方法和基于外部DEM數據的方法。第一類方法中，文獻[5]在雙天線交軌干涉中對主、輔圖像分別劃分兩個距離子帶構成的差分干涉對來計算絕對相位。為了利用同一場景不同干涉數據之間的冗余性，文獻[6]利用機載對飛模型的重疊區域建立的聯合相位偏移函數來估計絕對相位；文獻[7]通過聯合多基線干涉來求解相位偏移量。第一類方法需要滿足同一場景數據的冗余性，這限制了算法的普適性。第二類方法中，Perna針對機載雙天線干涉SAR提出了利用外源DEM進行兩步估計的方法[8]，首先利用粗精度DEM中提取的控制點估計初步的相位偏移量，然后利用初估計出的相位偏移量反演得到精度較低的高程數據，再利用控制點處的高程差與相位滿足的線性方程，采用最小二乘算法估計出第一步的估計誤差，通過第二步的迭代估計從而提高估計精度。

然而，重軌干涉SAR由于兩次飛行相互獨立，因殘余運動誤差無法完全抵消而引入未知的時變基線誤差。時變基線誤差會改變算法[8]第二步中建立的模型與線性求解的匹配性能，從而消弱算法估計精度甚至導致算法的不適用。受限于目前慣性導航和差分GPS精度的限制，時變基線誤差只能從數據中去估計。本文首先分析了時變基線誤差對利用外源DEM估計相位偏移量的影響，然后提出嵌入了時變基線誤差估計和補償的干涉相位偏置估計算法，從而改善了算法在重軌干涉中的適用性，最后通過機載C波段重軌干涉DEM反演實驗，驗證了算法的有效性。

1 機載干涉SAR相位偏移量估計模型

令相位解纏后的相位為φuw， 與絕對相位φ之間存在相位偏置φoff：

φoff=φuw-φ

(1)

利用地面控制點(Ground Control Point, GCP)去求解時，需獲得GCP較準確的地理信息，一般采用人為部署的角反射器。考慮從外源DEM提取的控制點作為定標點，外源DEM通過反向編碼投影到機載斜地幾何下，通過適當插值處理可以得到整個實驗場景的高程信息，然后提取一定數量的控制點用于計算相位偏置量。設控制點數量為N(N>1)，利用幾何關系可以分別計算出控制點的絕對相位值φGCP，式(1)可以寫成向量形式，即

φuw-1φoff=φGCP

(2)

(3)

式中，N×N矩陣W=diag(m1,…,mk,…,mN)，mk∈{0,1}表示掩膜標志。當控制點落入無效的掩膜區域時，其值mk=0，表明該點對估計沒有貢獻。相位偏置估計誤差為

(4)

圖1 利用外源DEM估計干涉相位偏置的流程圖

由于重軌干涉SAR高程靈敏度方程中，Δh與β(x,r)具有如下線性關系：

Δh=β(x,r)·Δφ(x,r)

(5)

式中：

(6)

式(5)反映了地形高程變化所引起的相位變化存在一定的比例關系。λ為波長，r為斜距，B為基線，B⊥為垂直基線，θ為下視角。

基于式(5)和式(6)，利用這個斜坡相位關系得到

(7)

Θ(x,r)≈h(x,r,φoff)-hGCP(x,r)

(8)

將式(7)寫成向量的形式，表示為

Δ=β·ε+Θ

(9)

由式(8)可知Θ是一個非零均值的矢量，此處將Θ寫成均值項v和非零的矢量Θ0，即

Θ=1v+Θ0

(10)

式中，v表示利用InSAR重建的DEM與外源DEM的相對偏量。將式(10)代入式(9)可得

Δ=β·ε+1v+Θ0

(11)

所以，式(11)可以表達為一個線性模型：

(12)

(13)

其加權最小二乘解為

(14)

式中，N×N矩陣W=diag(m1,…,mk,…,mN)，mk∈{0,1}表示掩膜標志。利用PBE和StopBE兩步估計便可以估計出最終的相位偏置。

值得注意的是，式(12)能夠采用最小二乘求解的前提是β-Δ具備線性性。實際上，β-Δ擬合的直線非常微弱地依賴于外源DEM的隨機波動誤差[8]。導致外源DEM波動誤差與β是相互獨立的，特別是當利用機載InSAR產生的DEM與外源DEM彼此獨立的情況。這并不意味著向量β的方差越小， StopBE估計就越準確。從式(13)可以看出，估計誤差是將式(12)中的噪聲Θ0投影到由系統矩陣H的列向量所展開的子空間，這個投影與HTH的行列式呈反比例關系：

(15)

式(15)表明估計誤差與β的方差呈比例關系，其比例常數為N2。更重要的是，式(15)中β的方差越大，式(13)引入的估計誤差就越小。因此， StopBE特別適合于機載干涉，因為機載干涉SAR具有較寬的下視角θ，β方差就特別大。

2考慮時變基線誤差補償的相位偏置量估計模型

2.1 時變基線誤差對相位偏置量估計算法的影響

時變基線是完全隨機的，本小節進行數值仿真實驗。雷達參數和基線信息如表1所示，仿真場景區的高程如圖2(a)所示，依據雷達參數和場景高程得到的絕對干涉相位如圖2(b)所示。仿真試驗時在圖2(b)中人為加入相位偏置量。由于外源DEM存在誤差，從DEM提取出的控制點的高程誤差可視為隨機噪聲，假設其滿足均值為μ、方差為σ2的高斯隨機分布G(μ,σ2)，并按照均勻分布原則選取不同位置的控制點。通過控制μ和σ2，從而模擬不同精度的DEM誤差。

表1 仿真時雷達參數和基線信息

(a)場景高程

(b)根據雷達參數和地形仿真的絕對相位圖2仿真場景的高程圖及干涉絕對相位

為了說明時變基線的影響，實驗中增加了時變基線誤差。所加入的時變基線誤差的水平分量dx和豎直分量dz如圖3(a)所示。時變基線誤差引起的每個像素的相位誤差如圖3(b)所示，可以看到時變基線誤差表現為相位波動。

為了定量分析時變基線對PBE和StopBE估計的影響，分別仿真了不同精度的DEM下的情況，設置了μ=10 m， σ=0, 0.5,…,5 m的情況和σ=2 m， μ=-10,-9,…,10 m的情況，從DEM中提取了1 000個地面控制點用于估計相位偏置，并給出了利用PBE和StopBE估計相位偏置反演出的高程誤差隨σ和μ的變化關系。圖4為μ=10 m，σ=0, 0.5,…,5 m的實驗結果，其中圖4(a)為沒有添加時變基線誤差的結果，圖4(b)為添加了時變基線誤差的結果。圖5為σ=2 m，μ=-10,-9,…,10 m的實驗結果，圖5(a)和圖5(b)分別對應沒有添加時變基線誤差和添加了時變基線誤差下的估計結果。從實驗結果可以看出，盡管DEM噪聲對PBE估計影響很大，但在不存在時變基線誤差的情況下，選擇一定數量的GCP點，通過PBE+StopBE的迭代估計的結果總能達到一個比較小的誤差；然而，在存在時變基線誤差的情況下，PBE+StopBE算法無法獲得一個較好的結果。因此，相位偏置估計算法在重軌干涉中運用時必須考慮時變基線誤差的補償。

(a)添加的時變基線誤差

(b)時變基線誤差導致的相位誤差圖3 仿真添加的時變基線誤差分量及其導致的相位誤差

(a)不存在時變基線誤差的情況

(b)存在時變基線誤差的情況圖4μ=10 m, σ=0, 0.5,…,5 m時PBE和StopBE估計相位偏置反演出的高程誤差隨σ的變化關系

(a)不存在時變基線誤差的情況

(b)存在時變基線誤差的情況圖5μ=-10,-9,…,10 m, σ=2 m時PBE和StopBE估計相位偏置反演出的高程誤差隨μ的變化關系

本文第1節分析過，StopBE采用最小二乘估計的前提是β-Δ的線性性能夠得到保持。圖6是不存在時變基線誤差(圖6(a))和存在時變基線誤差(圖6(b))時的β-Δ的關系及其擬合的直線。可以看出，不存在時變基線誤差的情況下，β-Δ的線性性能夠得到很好的保持。當時變基線誤差存在時β-Δ的線性性變差，從而影響了StopBE估計結果。由于時變基線是完全隨機和未知的，此處無法寫出其具體的解析表達式。然而，可以通過從數據域中估計并補償時變基線誤差，減小其對β-Δ的影響，以確保相位偏置量估計算法能夠使用。

(a)不存在時變基線誤差的情況

(b)存在時變基線誤差的情況

圖6不存在和存在時變基線誤差時的β-Δ的關系

2.2 時變基線誤差估計與補償算法

多斜視(Multisquint)方法[9]估計時變基線誤差時利用了方位配準誤差與不同子視圖像差分相位的關系。為了繞開低相干區目標對時變基線誤差估計的影響，改進的多斜視(Refined Multisquint)方法[10]被提出并成功用于估計時變基線誤差。對經過配準后的主S1、輔S2圖像分別劃分K個子孔徑圖像，通過計算相鄰子孔徑的干涉相位差，可以得到K-1個差分干涉相位。第k個差分干涉相位表達為

(16)

對差分干涉相位進行加權平均就得到時變基線的一階導數[10]：

(18)

Abxz=blos

(19)

式中：

(20)

(21)

(22)

其加權最小二乘解為

bxz=(ATWA)-1ATWblos

(23)

式中，加權矩陣W為

(24)

圖7 時變基線的距離空變幾何模型

通過對計算出的時變基線導數進行積分可以求解出時變基線誤差，未知的時變基線常數項可以利用外部DEM來消除，然后把計算出的時變基線誤差補償到最終的干涉相位上。

2.3 嵌入時變基線誤差補償的相位偏置量估計算法

重軌干涉SAR系統雙通道重復飛行的運動誤差相互獨立，無法完全相互抵消，因而引入了未知的時變基線誤差，所以必須采用殘余運動誤差補償算法，減小時變基線誤差φΔ對相位偏置估計算法的影響。圖8是嵌入了時變基線誤差估計的相位偏置估計算法。首先通過多斜視算法估計時變基線誤差，并對相位進行修正，然后用修正之后的相位結果進行相位偏置估計。該算法也可以迭代進行，以提高估計精度。

圖8 考慮時變基線誤差估計的相位偏置估計算法流程(虛線框為時變基線估計算法)

時變基線導數對頻譜間隔Δfsub的敏感度因子為

(25)

時變基線誤差的導數估計的標準差為

(26)

由此可知： 1) 劃分子視圖像數量K越多，越有利于提高時變基線導數估計的精度，然而劃分數量過多，會導致圖像信噪比變低，相干系數變低，從而使得差分干涉相位噪聲標準差增大，故實際處理時要綜合考慮； 2) 頻譜間隔越大，越有利于提高時變基線導數的估計精度，但頻譜間隔過大，會導致劃分的子孔徑帶寬偏小，會使干涉相位圖中的噪聲標準差變大，需要權衡考慮。

3 實測數據處理

為了驗證本文算法的有效性，對中國科學院電子學研究所的C波段機載重軌干涉SAR系統獲取的重軌干涉SAR數據進行處理。重軌干涉SAR參數如表2所示，成像算法采用了Ewk+DMA[11]+PTA[12]算法以減小運動補償帶來的誤差量，得到了聚焦良好的SAR圖像。

表2 重軌干涉SAR參數

處理場景中包括了3個定標點(A,B和C)，場景幅度圖和定標點位置如圖9所示。實驗中使用的DEM數據為垂直精度20 m、水平精度30 m的Aster DEM數據。Aster DEM轉化為對應場景高程數據，如圖10所示。圖像的多普勒有效帶寬為270 Hz，為了兼顧圖像信噪比和時變基線導數估計的精度，對配準之后的主、輔圖像分別劃分7個方位向子孔徑，每個子孔徑為38.57 Hz。采用本文方法估計和補償時變基線誤差，估計出來的時變基線的水平和豎直分量如圖11所示，估計出的殘余運動基本上在mm量級。

時變基線誤差補償之后，干涉處理時對方位進行了四視處理，去除平地并相位濾波之后的干涉結果如圖12所示。為了比較方法的有效性，此處對比了3種算法的結果：第一，先利用3個角反射器A,B,C計算出了相位偏置，然后反演高程，結果如圖13所示；第二，利用Aster DEM提取均勻分布的2 000個GCP點，并利用PBE+StopBE算法，迭代3次計算出相位偏置；第三，利用考慮了時變基線誤差估計的相位偏置估計方法，迭代3次計算出相位偏置，最后反演出的高程如圖14所示。3種情況估計出的相位偏置及在A,B,C定標點處反演的高程結果如表3所示。

通過比較實驗結果可以發現：直接利用定標點估計相位偏置，反演出的高程精度在2 m左右；直接利用Aster DEM提取控制點，采用PBE+StopBE迭代估計相位偏置，反演出的高程精度在6 m以內；利用本文改進的方法，即考慮了時變基線誤差估計和補償的相位偏置估計方法，最終反演出的高程誤差在4 m以內。改進后的方法減小了時變基線誤差對無人工控制點相位偏移量估計算法的影響，提高了DEM重建精度，適用于場景中無定標點的機載InSAR的DEM反演。

圖9 實驗場景幅度圖(A,B,C分別表示3個角反射器定標點)

圖10 從Aster DEM 轉化過來的對應成像場景的高程圖

圖11 估計到的時變基線水平和豎直分量

圖12 去平地后的干涉相位

圖13 利用角反射器估計相位偏置最后反演出來的高程圖

圖14 利用Aster-DEM數據按照本文方法計算相位偏置并反演的高程圖

表3 實驗結果比較

4 結束語

本文詳細分析了時變基線誤差對機載重軌干涉SAR相位偏置量估計算法的影響，指出了時變基線誤差會降低估計模型與線性求解的匹配性能，從而消弱算法估計精度的問題。針對這個問題，本文提出了嵌入時變基線誤差估計和補償的相位偏置量估計算法。最后利用機載C波段0.5 m高分辨率重軌干涉SAR實測數據驗證了該算法的有效性，并結合地面定標點驗證了該算法的高程重建精度在4 m以內。該方法簡單快速，消除了機載重軌干涉SAR獲取DEM時對人工定標點的依賴性，適用于無定標點情況下機載InSAR的DEM反演。

[1] 徐華平,陳杰,周蔭清,等. 干涉SAR三維地形成像數據處理技術綜述[J]. 雷達科學與技術, 2006, 4(1):15-21.

XU Huaping, CHEN Jie, ZHOU Yinqing, et al. A Survey of Interferometric SAR Topography Mapping Data Processing Technique[J]. Radar Science and Technology, 2006, 4(1):15-21.(in Chinese)

[2] 方東生,呂孝雷,李緣廷,等. 運動補償對機載SAR重軌干涉成像的影響分析[J]. 雷達科學與技術, 2016, 14(4):355-363.

FANG Dongsheng, LYU Xiaolei, LI Yuanting, et al. Effect of Motion Compensation on Airborne Repeat Pass InSAR Imaging[J]. Radar Science and Technology, 2016, 14(4):355-363.(in Chinese)

[3] NEUMANN M, HENSLEY S, AHMED R, et al. UAVSAR PolInSAR and Tomographic Products for Natural Media Characterization[C]∥10th European Conference on Synthetic Aperture Radar, Berlin, Germany:VDE, 2014:225-228.

[4] TEBALDINI S, NAGLER T, ROTT H, et al. Imaging the Internal Structure of an Alpine Glacier via L-Band Airborne SAR Tomography [J]. IEEE Trans on Geoscience and Remote Sensing, 2016, 54(12):7197-7209.

[5] MADSEN S N, ZEBKER H A. Automated Absolute Phase Retrieval in Across-Track Interferometry[J]. International Geoscience and Remote Sensing Sympo-sium, Houston, TX:IEEE, 1992:1582-1584.

[6] MURA J C, PINHEIRO M, ROSA R, et al. A Phase-Offset Estimation Method for InSAR DEM Generation Based on Phase-Offset Functions[J]. Remote Sensing, 2012, 4(3):745-761.

[7] GATTI G, TEBALDINI S, D’ALESSANDRO M M, et al. ALGAE:A Fast Algebraic Estimation of Interferogram Phase Offsets in Space-Varying Geometries [J]. IEEE Trans on Geoscience and Remote Sensing, 2011, 49(6):2343-2353.

[8] PERNA S, ESPOSITO C, BERARDINO P, et.al. Phase Offset Calculation for Airborne InSAR DEM Generation Without Corner Reflectors[J]. IEEE Trans on Geoscience and Remote Sensing, 2015, 53(5):2713-2726.

[9] SCHEIBER R, MOREIRA A. Coregistration of Interferometric SAR Images Using Spectral Diversity[J]. IEEE Trans on Geoscience and Remote Sensing, 2000, 38(5):2179-2191.

[10] REIGHER A, PRATS P, MALLORQUI J J. Refined Estimation of Time-Varying Baseline Errors in Airborne SAR Interferometry [J]. IEEE Geoscience and Remote Sensing Letters, 2006,3(1):145-149.

[11] MENG D, HU D, DING C. A New Approach to Airborne High Resolution SAR Motion Compensation for Large Trajectory Deviations [J]. Chinese Journal of Electronics, 2012, 21(4):764-769.

[12] DE MACEDO K A C, SCHEIBER R. Precise Topography and Aperture Dependent Motion Compensation for Airborne SAR[J]. IEEE Geoscience and Remote Sensing Letters, 2005, 2(2):172-176.