地形高程自適應的星載InSAR圖像配準方法

張金強， 索志勇， 李真芳， 保　錚

(西安電子科技大學 雷達信號處理國家重點實驗室，陜西 西安 710071)

星載干涉合成孔徑雷達(Interferometric Synthetic Aperture Radar， InSAR)技術利用不同觀測幾何下，獲取的兩幅合成孔徑雷達(Synthetic Aperture Radar，SAR)復圖像間的干涉相位反演來獲取地表三維地形[1-3]．為了從SAR復圖像對中獲取干涉相位信息，需要對其進行配準處理，以確保用于干涉處理的圖像像素對應于地面同一散射單元[4-5]．圖像配準精度將影響InSAR系統的最終產品精度[6]．因此，研究穩健高精度InSAR圖像配準方法具有重要意義．已有的InSAR圖像配準方法主要可分為兩大類：第1類方法首先利用基于圖像數據[7-8]或圖像特征[9-10]的方法估計控制點處配準偏移量，然后根據二維二元低階多項式模型計算其他像素處的配準偏移量．該類方法沒有考慮地形高程對配準偏移量的影響，當基線較長或觀測場景內地形復雜度較高時，利用二維二元低階多項式模型擬合配準偏移量的精度較低[11]；當利用控制點處配準偏移量估計多項式系數時，控制點數量和位置分布對系數估計精度的影響較大；第2類方法首先利用幾何配準法[12]計算每一像素的初始配準偏移量，然后利用控制點處高精度配準偏移量根據二維二元一階多項式模型校正初始配準偏移量[11]，從而獲取全場景高精度配準偏移量．該類方法需要進行逐像素幾何配準處理，運算效率較低．

針對上述問題，筆者提出了地形高程自適應的降維圖像配準方法．首先，給出了主輔圖像間的配準偏移量計算模型；然后，為了追蹤配準偏移量隨地形高程的變化，根據配準偏移量計算模型提出了包含地形高程項的二維三元一階多項式模型．利用該模型不僅可以精確擬合主輔圖像間的配準偏移量，而且能夠避免逐像素幾何配準處理；最后，通過降維處理將需要利用控制點處配準偏移量估計的多項式系數減少為兩個常數項系數，提高了系數估計的穩健性．實測數據處理結果驗證了文中方法的精確性和穩健性．

1　星載InSAR降維圖像配準方法

1.1　配準偏移量計算模型

主輔圖像間的配準偏移量計算流程包括兩步:首先，通過正向定位獲取主圖像像素對應的目標三維位置;然后，通過反向定位獲取上述目標在輔圖像上的坐標．針對主圖像上某一像素(am，rm)，將成像參數和軌道數據代入下列3式：

所示的距離-多普勒模型[11-13]求解目標三維位置，即正向定位，其中，R表示雷達斜距；t表示方位時間；P(t)= (Px，Py，Pz)T，表示衛星位置矢量；T= (Tx，Ty，Tz)T，表示目標三維位置；fdi表示成像多普勒頻率；V(t)= (Vx，Vy，Vz)T，表示衛星速度矢量；λ表示雷達波長；h表示目標高程；P(t)和V(t)隨時間變化；Re和Rp分別表示赤道和極地半徑，且Rp= (1-f)(Re+h)，f為參考橢球扁率；P、V和T均定義在世界測地系統(World Geodetic System-1984，WGS-84)坐標系下．

利用目標三維位置T、輔圖像成像參數和軌道數據，首先根據式(2)計算該目標對應的輔圖像方位時間t，然后根據式(1)計算時刻t目標到雷達的斜距R，最后根據

計算該目標在輔圖像上的二維坐標(as，rs)，其中，t0表示圖像方位向第1個像素對應的方位時間，fp表示脈沖重復頻率，R0表示圖像距離向第1個像素對應的雷達斜距，fs表示脈沖采樣頻率，c表示光速．上述過程稱為反向定位．

1.2　配準偏移量擬合模型

由配準偏移量計算模型可知，目標在輔圖像上的方位和距離坐標隨其在主圖像上的方位、距離坐標和高程變化．傳統配準方法假設目標高程對配準偏移量的影響可以忽略，用含有方位和距離坐標項的二維二元低階多項式模型擬合配準偏移量隨主圖像方位和距離坐標的變化．當主輔圖像間的基線較長或觀測場景內地形復雜度較高時，地形高程對配準偏移量的影響不能忽略．為了追蹤配準偏移量隨地形高程的變化，可用

所示的包含地形高程項的二維三元一階多項式模型擬合主輔圖像間的坐標關系，其中，d0～d3和g0～g3表示擬合系數．

1.3　系數估計

式(6)和式(7)所示的二維三元一階多項式模型存在8個未知系數．傳統配準方法利用控制點處配準偏移量估計上述未知系數，控制點數量和位置分布對系數估計精度影響較大．為了降低控制點對系數估計的影響，提高系數估計的穩健性，文中給出了利用主輔圖像的成像參數和軌道數據計算一次項系數的公式．一次項系數計算公式的推導包括三步:首先，根據正向定位模型，推導目標三維位置相對于主圖像方位時間、雷達斜距和目標高程的一階偏導數;然后，根據反向定位模型，推導目標對應輔圖像方位時間和雷達斜距相對于目標三維位置的一階偏導數;最后，結合上述兩步結果得到一次項系數的計算公式．

根據正向定位模型，將式(1)～式(3)中目標三維位置分別對方位時間、雷達斜距和目標高程求一階偏導數，得到

其中，Am表示主衛星加速度矢量．其他參數與1.1節的參數含義相同，下標m表示主圖像對應的參數．

根據反向定位模型，將式(1)和式(2)中方位時間和雷達斜距分別對目標三維位置求一階偏導數，得到

其中，As表示輔衛星加速度矢量．其他參數與1.1節的參數含義相同，下標s表示輔圖像對應的參數．

根據式(4)～式(5)和式(8)～式(14)，得到目標在輔圖像上的二維坐標相對于其在主圖像上的二維坐標和高程的一階偏導數為

其中，d3和g3的單位為 像素/m．上述偏導數隨主圖像二維坐標和目標高程變化．由于星載SAR的雷達斜距較大[14]，在一定觀測場景范圍內，上述偏導數的變化對配準偏移量的影響可以忽略．因此，在實測數據處理過程中，利用觀測場景中心處目標對應的主輔圖像的成像參數和軌道數據，根據式(15)～式(20)計算d1～d3和g1～g3．

圖1 文中方法的處理流程

1.4　處理流程

如圖1所示，地形高程自適應的降維圖像配準方法處理流程包括如下步驟：

(1) 利用觀測場景中心處目標對應的主輔圖像的成像參數和軌道數據，根據式(15)～式(20)計算d1～d3和g1～g3．

(2) 在主圖像上選取一定數量的控制點，采用基于圖像數據或圖像特征的方法估計控制點處配準偏移量．

(3) 根據式(6)和式(7)，補償控制點處配準偏移量的一次分量，利用剩余配準偏移量估計d0和g0．

(4) 根據式(6)和式(7)，計算其他像素處配準偏移量．

文中方法需要輸入主圖像像素對應的目標高程，其可以利用先驗數字高程模型(Digital Elevation Model，DEM)通過正向定位獲取．

2　實測數據處理與分析

利用TerraSAR-X分別于2010年3月8日和2010年3月19日對某地區進行重復航過觀測任務時，獲取的兩幅SAR復圖像來驗證文中方法，取2010年3月8日獲取的SAR復圖像為主圖像．主要系統參數如表1所示，SAR幅度圖像如圖2(a)所示．在圖2(a)方框所示的區域內選擇12個特顯點以定量評估配準精度，如圖2(b)所示．圖2(c)為觀測場景對應的航天飛機雷達地形測繪任務(Shuttle Radar Topography Mission，SRTM)DEM，已將其投影到主圖像斜距平面，觀測場景幅寬為 15 km× 15 km，觀測場景內目標高程變化為 695 m．

表1　TerraSAR-X主要系統參數

圖2　TerraSAR-X重復航過數據

利用圖1所示的處理流程對主輔圖像進行配準處理，在主圖像內均勻選擇100個控制點．控制點處配準偏移量采用實相關函數法[7]估計，圖像子塊的大小為 64像素× 64像素．配準后主輔圖像間的干涉相位圖如圖3所示．利用文中方法配準后的主輔圖像可以獲得清晰的干涉相位圖，驗證了文中方法的有效性．

為了驗證文中方法的精確性和穩健性，利用蒙特卡羅實驗對比分析文中方法與二階多項式方法(利用二維二元二階多項式模型擬合配準偏移量，并且多項式系數均利用控制點處配準偏移量估計)在不同控制點數量下的配準性能．蒙特卡羅實驗次數為100次．在主圖像內均勻選擇900個控制點，控制點處配準偏移量采用實相關函數法估計，圖像子塊的大小為 64像素× 64像素．每次實驗從上述900個控制點中隨機選擇一定數量的控制點估計多項式擬合系數．配準誤差利用圖2(b)所示的12個特顯點進行估計．對配準后主輔圖像上的同名特顯點進行200倍插值，將峰值間的距離向和方位向偏移量作為該特顯點處的配準誤差．圖4給出了二階多項式方法和文中方法配準誤差均方根值隨控制點數量的變化．由圖4可知，二階多項式方法配準誤差隨控制點數量變化較大，而文中方法配準誤差隨控制點數量變化較小．例如，當控制點數量為10個時，二階多項式方法的方位向和距離向配準誤差分別為0.79像素和0.86像素，文中方法的方位向和距離向配準誤差分別為0.05像素和0.07像素; 當控制點數量為100個時，二階多項式方法的方位向和距離向配準誤差分別為0.05像素和0.12像素，文中方法的方位向和距離向配準誤差分別為0.04像素和0.05像素．由表1可知，主輔圖像間的垂直基線大于水平基線，導致目標高程對距離向配準偏移量的影響大于方位向配準偏移量，所以二階多項式方法的距離向配準誤差大于方位向配準誤差．文中方法可以追蹤配準偏移量隨目標高程的變化，所以文中方法的方位向和距離向配準誤差近似相同，并且可以滿足0.1像素的配準精度要求．上述結果驗證了文中方法的精確性和穩健性．

圖3 干涉相位圖圖4 二階多項式方法和文中方法配準誤差統計結果

為了直觀證明文中方法相對于二階多項式方法的精確性和穩健性，圖5給出了控制點數量為10個時，100次試驗的12個特顯點處的配準誤差分布．圖5(a)～(d)分別為二階多項式方法的方位向和距離向配準誤差、文中方法的方位向和距離向配準誤差．由圖5可知，二階多項式方法的方位向和距離向配準誤差最大值分別為3.47像素和3.45像素，而文中方法的方位向和距離向配準誤差最大值分別為0.17像素和0.30像素．

圖5　二階多項式方法和文中方法配準誤差分布

3　結 束 語

針對長基線或復雜場景星載InSAR圖像配準問題，文中提出了地形高程自適應的降維圖像配準方法．首先，給出了主輔圖像間的配準偏移量計算模型；然后，提出了包含地形高程項的二維三元一階多項式模型．該模型可以追蹤配準偏移量隨地形高程的變化，提高了配準偏移量擬合精度；最后，利用主輔圖像的成像參數和軌道數據計算多項式一次項系數，將需要利用控制點處配準偏移量估計的系數減少為兩個常數項系數，提高了系數估計的穩健性．實測數據處理結果驗證了文中方法的精確性和穩健性．

參考文獻：

[1] YAGUE-MARTINEZ N, PRATS-IRAOLA P, GONZALEZ F R, et al. Interferometric Processing of Sentinel-1 TOPS Data[J]. IEEE Transactions on Geoscience and Remote Sensing, 2016, 54(4): 2220-2234.

[2]LIANG C, FIELDING E J. Interferometry with ALOS-2 Full-aperture ScanSAR Data[J]. IEEE Transactions on Geoscience and Remote Sensing, 2017, 55(5): 2739-2750.

[3]于瀚雯, 保錚. 利用L1范數的多基線InSAR相位解纏繞技術[J]. 西安電子科技大學學報, 2013, 40(4): 37-41.

YU Hanwen, BAO Zheng.L1-norm Method for Multi-baseline InSAR Phase Unwarpping[J]. Journal of Xidian University, 2013, 40(4): 37-41.

[4]靳峰, 馮大正. 一種快速準確的圖像配準算法[J]. 西安電子科技大學學報, 2015, 42(6): 88-93.

JIN Feng, FENG Dazheng. Fast and Precise Image Registration Algorithm[J]. Journal of Xidian University, 2015, 42(6): 88-93.

[5]HAN Y, BOVOLO F, BRUZZONE L. Segmentation-based Fine Registration of Very High Resolution Multitemporal Images[J]. IEEE Transactions on Geoscience and Remote Sensing, 2017, 55(5): 2884-2897.

[6]KRIEGER G, MOREIRA A, FIEDLER H, et al. TanDEM-X: a Satellite Formation for High-resolution SAR Interferometry[J]. IEEE Transactions on Geoscience and Remote Sensing, 2007, 45(11): 3317-3341.

[7]STONE H S, ORCHARD M T, CHANG E C, et al. A Fast Direct Fourier-based Algorithm for Subpixel Registration of Images[J]. IEEE Transactions on Geoscience and Remote Sensing, 2001, 39(10): 2235-2243.

[8]薛海偉, 馮大正. 一種快速的InSAR圖像精配準方法[J]. 西安電子科技大學學報, 2016, 43(3): 172-178.

XUE Haiwei, FENG Dazheng. Fast Subpixel Registration Method for InSAR Images[J]. Journal of Xidian University, 2016, 43(3): 172-178.

[9]SUN K, LIU L M, TAO W B. Image Matching via Feature Fusion and Coherent Constraint[J]. IEEE Geoscience and Remote Sensing Letters, 2017, 14(3): 289-293.

[10]FAN J W, WU Y, WANG F, et al. New Point Matching Algorithm Using Sparse Representation of Image Patch Feature for SAR Image Registration[J]. IEEE Transactions on Geoscience and Remote Sensing, 2017, 55(3): 1498-1510.

[11]NITTI D O, HANSSEN R F, REFICE A, et al. Impact of DEM-assisted Coregistration on High-resolution SAR Interferometry[J]. IEEE Transactions on Geoscience and Remote Sensing, 2011, 49(3): 1127-1143.

[12]SANSOSTI E, BERARDINO P, MANUNTA M, et al. Geometrical SAR Image Registration[J]. IEEE Transactions on Geoscience and Remote Sensing, 2006, 44(10): 2861-2870.

[13]劉艷陽, 李真芳, 楊娟娟, 等. 分布式衛星InSAR目標定位近似閉式解[J]. 西安電子科技大學學報, 2012, 39(4): 87-93.

LIU Yanyang, LI Zhenfang, YANG Juanjuan, et al. Quasi-closed-form Solution for Distributed Satellites InSAR Geolocation[J]. Journal of Xidian University, 2012, 39(4): 87-93.

[14]李錦偉, 李真芳, 侯英龍, 等. 基于遞推公式的星載SAR高效正向定位算法[J]. 電子與信息學報, 2014, 36(2): 409-414.

LI Jinwei, LI Zhenfang, HOU Yinglong, et al. A Novel Efficient Spaceborne SAR Geolocation Method Based on Recursion Formulae[J]. Journal of Electronics & Information Technology, 2014, 36(2): 409-414.