張金強(qiáng)， 索志勇， 李真芳， ?！″P

(西安電子科技大學(xué) 雷達(dá)信號(hào)處理國家重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710071)

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

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

1　星載InSAR降維圖像配準(zhǔn)方法

1.1　配準(zhǔn)偏移量計(jì)算模型

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

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

利用目標(biāo)三維位置T、輔圖像成像參數(shù)和軌道數(shù)據(jù)，首先根據(jù)式(2)計(jì)算該目標(biāo)對(duì)應(yīng)的輔圖像方位時(shí)間t，然后根據(jù)式(1)計(jì)算時(shí)刻t目標(biāo)到雷達(dá)的斜距R，最后根據(jù)

計(jì)算該目標(biāo)在輔圖像上的二維坐標(biāo)(as，rs)，其中，t0表示圖像方位向第1個(gè)像素對(duì)應(yīng)的方位時(shí)間，fp表示脈沖重復(fù)頻率，R0表示圖像距離向第1個(gè)像素對(duì)應(yīng)的雷達(dá)斜距，fs表示脈沖采樣頻率，c表示光速．上述過程稱為反向定位．

1.2　配準(zhǔn)偏移量擬合模型

由配準(zhǔn)偏移量計(jì)算模型可知，目標(biāo)在輔圖像上的方位和距離坐標(biāo)隨其在主圖像上的方位、距離坐標(biāo)和高程變化．傳統(tǒng)配準(zhǔn)方法假設(shè)目標(biāo)高程對(duì)配準(zhǔn)偏移量的影響可以忽略，用含有方位和距離坐標(biāo)項(xiàng)的二維二元低階多項(xiàng)式模型擬合配準(zhǔn)偏移量隨主圖像方位和距離坐標(biāo)的變化．當(dāng)主輔圖像間的基線較長或觀測場景內(nèi)地形復(fù)雜度較高時(shí)，地形高程對(duì)配準(zhǔn)偏移量的影響不能忽略．為了追蹤配準(zhǔn)偏移量隨地形高程的變化，可用

所示的包含地形高程項(xiàng)的二維三元一階多項(xiàng)式模型擬合主輔圖像間的坐標(biāo)關(guān)系，其中，d0～d3和g0～g3表示擬合系數(shù)．

1.3　系數(shù)估計(jì)

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

根據(jù)正向定位模型，將式(1)～式(3)中目標(biāo)三維位置分別對(duì)方位時(shí)間、雷達(dá)斜距和目標(biāo)高程求一階偏導(dǎo)數(shù)，得到

其中，Am表示主衛(wèi)星加速度矢量．其他參數(shù)與1.1節(jié)的參數(shù)含義相同，下標(biāo)m表示主圖像對(duì)應(yīng)的參數(shù)．

根據(jù)反向定位模型，將式(1)和式(2)中方位時(shí)間和雷達(dá)斜距分別對(duì)目標(biāo)三維位置求一階偏導(dǎo)數(shù)，得到

其中，As表示輔衛(wèi)星加速度矢量．其他參數(shù)與1.1節(jié)的參數(shù)含義相同，下標(biāo)s表示輔圖像對(duì)應(yīng)的參數(shù)．

根據(jù)式(4)～式(5)和式(8)～式(14)，得到目標(biāo)在輔圖像上的二維坐標(biāo)相對(duì)于其在主圖像上的二維坐標(biāo)和高程的一階偏導(dǎo)數(shù)為

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

圖1 文中方法的處理流程

1.4　處理流程

如圖1所示，地形高程自適應(yīng)的降維圖像配準(zhǔn)方法處理流程包括如下步驟：

(1) 利用觀測場景中心處目標(biāo)對(duì)應(yīng)的主輔圖像的成像參數(shù)和軌道數(shù)據(jù)，根據(jù)式(15)～式(20)計(jì)算d1～d3和g1～g3．

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

(3) 根據(jù)式(6)和式(7)，補(bǔ)償控制點(diǎn)處配準(zhǔn)偏移量的一次分量，利用剩余配準(zhǔn)偏移量估計(jì)d0和g0．

(4) 根據(jù)式(6)和式(7)，計(jì)算其他像素處配準(zhǔn)偏移量．

文中方法需要輸入主圖像像素對(duì)應(yīng)的目標(biāo)高程，其可以利用先驗(yàn)數(shù)字高程模型(Digital Elevation Model，DEM)通過正向定位獲取．

2　實(shí)測數(shù)據(jù)處理與分析

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

表1　TerraSAR-X主要系統(tǒng)參數(shù)

圖2　TerraSAR-X重復(fù)航過數(shù)據(jù)

利用圖1所示的處理流程對(duì)主輔圖像進(jìn)行配準(zhǔn)處理，在主圖像內(nèi)均勻選擇100個(gè)控制點(diǎn)．控制點(diǎn)處配準(zhǔn)偏移量采用實(shí)相關(guān)函數(shù)法[7]估計(jì)，圖像子塊的大小為 64像素× 64像素．配準(zhǔn)后主輔圖像間的干涉相位圖如圖3所示．利用文中方法配準(zhǔn)后的主輔圖像可以獲得清晰的干涉相位圖，驗(yàn)證了文中方法的有效性．

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

圖3 干涉相位圖圖4 二階多項(xiàng)式方法和文中方法配準(zhǔn)誤差統(tǒng)計(jì)結(jié)果

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

圖5　二階多項(xiàng)式方法和文中方法配準(zhǔn)誤差分布

3　結(jié) 束 語

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

參考文獻(xiàn)：

[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范數(shù)的多基線InSAR相位解纏繞技術(shù)[J]. 西安電子科技大學(xué)學(xué)報(bào), 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]靳峰, 馮大正. 一種快速準(zhǔn)確的圖像配準(zhǔn)算法[J]. 西安電子科技大學(xué)學(xué)報(bào), 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圖像精配準(zhǔn)方法[J]. 西安電子科技大學(xué)學(xué)報(bào), 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]劉艷陽, 李真芳, 楊娟娟, 等. 分布式衛(wèi)星InSAR目標(biāo)定位近似閉式解[J]. 西安電子科技大學(xué)學(xué)報(bào), 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]. 電子與信息學(xué)報(bào), 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.